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ABSTRACT 

When  microbubbles  are  injected  in  a  turbulent  boundary  layer, 
they  promote  a  reduction  of  the  local  skin  friction  and  a  sharp 
mean  velocity  gradient  between  the  two  microbubble  free  regions 
in  the  boundary  layer,  i.e.,  the  outer  edge  and  the  near  wall 
region. 

In  connection  with  the  apparent  stability  of  the  bubbles  in 
the  boundary  layer,  the  stability  of  a  hollow  vortex  sheet  in  an 
infinite  expanse  of  fluid  as  well  as  in  the  presence  of  a  wall  are 
studied.  We  concluded  that  the  simple  hollow  vortex  sheet  is 
unconditionally  unstable  while  a  stability  range  was  found  for  a 
hollow  vortex  sheet  perturbed  In  the  presence  of  a  wall.  In  all 
cases,  however,  the  computations  suggest  that  the  hollow  vortex 
radius  is  invariably  a  destabilizing  factor. 

Regarding  the  second  asepct  of  the  problem,  a  heuristic, 
but  physically  sound  explanation  of  the  mechanisms  of  the  drag 
reduction  is  proposed.  Based  on  these  ideas  and  using  a 
phenomenlogical  approach,  the  flat  plate  boundary  layer  is 
numerically  analyzed. 
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CHAPTER  I 
INTRODUCTION 

1 . 1  Previous  Experimental  Work 

,  The  idea  of  reducing  the  skin  friction  between  water  and  a 
solid  by  injecting  air  close  to  the  surface  goes  back  to  Froude, 
1875.  It  was  hoped  that  the  lower  viscosity  of  the  in-jected 
fluid  close  to  the  surface  would  lead  to  a  reduction  in  the  skin 
friction.  Indeed,  from  the  very  beginning  experiments  validated 
this  idea.  Even  though  the  mechanisms  by  which  the  drag 
reduction  is  brought  about  are  not  yet  fully  understood,  the 
phenomenon  and  its  potential  applications  have  lured  the 
attention  of  physicists  and  engineers  ever  since. 

What  precludes  our  full  understanding  of  the  problem  is  that 
it  occurs  in  a  strictly  turbulent  flow.  Turbulence  itself  is 
included  in  the  category  of  one  of  the  major  unsolved  problems  of 
physics  today.  With  the  addition  of  small  bubbles  in  the  boundary 
layer,  the  drag  reduction  is  intimately  related  to  subtle  flow 
interactions  between  the  bubbles  and  the  carrier  fluid  in  the 
viscous  and  the  buffer  layers  of  the  boundary  layer.  Even  in  the 
case  of  the  simple  turbulent  boundary  layer,  today's  knowledge  of 
the  near  wall  region  is  far  from  complete  mainly  because  of 
experimental  difficulties  in  obtaining  reliable  data  close  to  the 
surface.  In  the  case  of  the  turbulent  boundary  layer  with 
microbubbles,  this  difficulty  is  further  aggravated  because  the 
thickness  of  the  film  of  fluid  between  the  bubble  layer  and  the 
wall  is  of  the  order  of  micrometers  and  the  bubbles  are  optically 
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The  interest  in  studying  this  phenomenon  '-as  intensified  in  the 
past  twenty  years  mainly  as  a  result  of  a  better  physical  under¬ 
standing  of  turbulent  flows  and  of  new  experimental  techniaues 
developed  to  study  more  complex  flowfields. 

Some  researchers  [1]  chose  to  generate  the  microbubbles  in  the 
boundary  layer  via  electrolysis.  This  approach,  however,  does  not 
seem  to  be  the  most  appropriate  if  one  is  interested  in  studying  the 
interactions  between  the  bubbles  and  the  mean  flow  in  the  boundary 
layer  because  of  the  presence  there  of  the  gas  generating  wire. 

The  study  of  the  turbulent  boundary  layer  flow  with  microbubble 
injection  was  undertaken  in  a  systematic  fashion  for  the  first  time 
in  the  Soviet  Union  [2,3,4,51  at  the  beginning  of  the  last  decade. 
These  researchers,  as  well  as  the  research  group  at  The  Pennsylvania 
State  University,  have  become  dedicated  more  recently  to  the  same 
problem  [6,7,8].  In  the  United  States,  this  group  chose  to  inject 
the  gaseous  fluid  in  the  boundary  layer  through  porous  plates. 

The  results  of  the  work  of  N.  K.  Madavan,  et  al.  [8]  as  well 
as  those  of  the  Soviet  researchers  show  that  the  local  skin 
friction  coefficient  can  be  reduced  to  within  20%  of  its 
undisturbed  value.  Figure  1.1  illustrates  some  of  these  findings. 

In  this  particular  work,  the  boundary  layer  was  observed  on  a 
flat  plate  at  zero  pressure  gradient.  In  the  absence  of  micro- 
bubbles,  the  boundary  layer  exhibited  the  classical  characteristics 
of  a  fully  developed  turbulent  boundary  layer,  as  is  evident  from 
the  data  which  are  reproduced  in  Figure  1.2  and  Figure  1.3. 
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As  reported  in  [8],  Figure  1.1  indicates  that  the  skin  friction 
is  decreased  as  soon  as  the  gas  starts  being  injected  through  the 
plate.  It  is  noted  also  that  the  amount  of  the  skin  friction 
reduction  generally  increases  with  increased  gas  flow  and  decreased 
free  stream  velocity.  Indeed,  the  results  suggest  that  for  a  given 
free  stream  velocity  there  is  an  optimum  gas  flow  rate  for  which  the 
skin  friction  reduction  is  a  maximum.  In  (8]  it  is  shown  that  the 
skin  friction  reduction  data  can  be  collapsed  on  the  same  curve  if 
plotted  against  the  ratio  of  the  gas  flow  to  the  total  flow  in  the 
boundary  layer. 

Generally,  for  any  gas  flow  rate,  it  was  observed  that  the 
maximum  skin  friction  reduction  occurred  immediately  downstream  of 
the  porous  section  and  that  the  skin  friction  gradually  recovers 
to  its  undisturbed  value  in  the  downstream  direction.  Madavan, 
et  al.  [6]  reports  that  the  skin  reduction  is  felt  as  far  down¬ 
stream  as  60  to  70  boundary  layer  thicknesses.  The  results  shown 
in  Figure  l.A  and  Figure  1.5,  taken  from  this  work,  are 
representative  of  the  conclusions  obtained. 

Both  N.  K.  Madavan,  et  al.  [6]  and  G.  S.  Migirenko,  et  al. 


[7]  report  that  at  the  higher  speed  flow  regimes  the  bubble  packing 
is  higher  while  the  bubble  mean  radius  is  smaller.  The  latter 


author  claims  that  the  highest  bubble  volumetric  concentrations 
of  the  order  of  80%  are  at  a  location  of  around  0.15*.  In  fact, 
contrary  to  the  findings  of  the  Soviet  researchers  in  [6],  it  is 
reported  that  the  bubble  radius  is  much  more  sensitive  to  the 
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flow  conditions  than  to  the  characteristics  of  the  porous  materials 
through  which  the  gas  is  injected.  This  seems  to  be  in  line  with 
experimental  results  reported  earlier  [9]. 

The  recovery  of  the  wall  mean  velocity  gradient  was  observed  to 
be  associated  with  the  migration  of  the  bubbles  away  from  the 
boundary  layer.  The  time  scale  of  the  migration  is  associated  with 
the  buoyant  force  and  is  consequently  of  longer  duration  for  the 
smaller  bubbles.  The  higher  packing  seems  to  be  directly  related 
to  the  extent  to  which  the  skin  friction  reduction  persists 
downstream. 

Another  interesting  feature  reported  in  both  [1]  and  [2]  is 
that  the  bubble  concentration  or,  in  other  words,  the  void  fraction 
dies  out  as  the  boundary  layer  edge  is  approached.  The  measurements 
of  the  mean  velocity  profile  in  this  part  of  the  boundary  layer 
show  that  it  is  not  affected  by  the  presence  of  the  bubbles  in  the 
boundary  layer.  The  LDA  measurements  in  this  region  of  the 
boundary  layer,  as  presented  in  (8),  are  reproduced  in  Figure  1.6. 

Finally,  both  the  American  and  Soviet  researchers  report  that 
the  high  frequency  signals  of  the  shear  stress  and  pressure 
fluctuations  close  to  the  surface  are  lost  with  the  injection  of 
the  gas. 

1 . 2  Problem  Statement 

The  problem  analyzed  in  the  following  chapters  has  two  distinct 
aspects,  i.e.,  the  stability  of  the  bubble  layer  in  the  shear  flow 
and  the  drag  reduction. 
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Figure  1.6.  Comparison  of  Mean  Velocicy  Profiles  in  the  Outer 
Part  of  the  Boundary  Layer  in  the  Presence  and 
Absence  of  Microbubbles. 


1.2.1  The  Stability 

In  light  of  the  brief  discussion  presented  above,  we  see  that 
the  insertion  of  the  bubble  layer  in  the  boundary  layer  appears  to 
give  rise  to  a  sharp  mean  velocity  variation  across  the  bubble 
layer.  The  above  statement  is  motivated  by  the  fact  that  while 
promoting  a  decrease  in  the  wall  velocity  gradient  the  bubbles  do 
not  affect  the  outer  region  of  the  boundary  layer. 

The  fact  that  the  bubble  layer  seems  to  be  stable  in  the 
turbulent  shear  flow  was  to  a  certain  extent  surprising  because  of 
Helmholtz's  classical  result  on  vortex  sheet  unconditional 
instability. 

The  question  of  the  stability  is  looked  upon  from  several 
points  of  view.  Given  the  two-dimensional  nature  of  the 
boundary  layer  flow,  the  stability  of  a  hollow  vortex  sheet  is 
considered  in  the  presence  and  in  the  absence  of  a  wall.  We 
study  the  problem  from  both  an  analytical  and  numerical  point 
of  view. 


1.2.2  The  Drag  Reduction 

As  is  clear  from  the  description  of  the  experimental  results, 
the  flowfield  is  rather  complex  and  the  mechanisms  by  which  the 
drag  reduction  is  brought  about  are  not  entirely  clear. 

Due  to  the  small  size  of  the  bubbles  and  the  two-dimensional 
nature  of  the  flow,  we  adopt  a  phenomenological  approach  in 
which  the  fluid  is  assumed  to  have  specially  varying  properties. 
The  obvious  limitations  of  this  model  are  discussed  in 
Chapter  IV. 


CHAPTER  II 


STABILITY  OF  A  HOLLOW  VORTEX  SHEET 

2 . 1  Introduction 

Helmholtz's  classical  result  on  the  unconditional  instability 
of  a  plane  velocity  discontinuity  has  been  re-examined  by  several 
authors  from  both  computational  [10,11,12)  and  analytical  [13] 
points  of  view. 

In  von  Karman's  classical  analysis  (see  Sir  H.  Lamb  [13]),  as 
well  as  in  all  subsequent  works,  the  plane  surface  velocity 
discontinuity  is  modeled  by  a  doubly  infinite  array  of  discrete 
vortices  the  circulation  of  which  depends  on  the  velocity  jump 
across  the  surface.  In  addition  to  corroborating  Helmholtz's 
result,  von  Karman,  in  his  elegant  analysis,  showed  also  that  the 
presence  of  a  wall  has  no  influence  on  the  vortex  sheet  so  far 
as  its  stability  is  concerned. 

In  the  following  sections,  we  try  to  devise  a  simple  model 
that  may  be  representative  of  the  two-phase  phenomena  observed  and 
will  give  us  a  working  base  for  the  studv  of  the  stnbilitv.  "Ihe 
two-dimensional  nature  of  the  boundary  layer  leads  one  naturally 
to  imagine  the  bubbles  as  cylinders.  In  the  realm  of  inviscid 
flows,  we  analyze  the  problem  of  the  stability  of  a  hollow  vortex 
sheet  from  two  different  points  of  view. 

In  the  following  sections,  we  take  the  spirit  of  von  Karman's 
work  and  extend  his  analysis  in  order  to  consider  hollow  core 
vortices  instead  of  point  vortices.  The  main  difference  becomes 
the  superposition  on  the  original  flowfield  of  the  appropriate 
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velocity  potentials  in  order  to  guarantee  the  existence  of 


circular  strearalinea  centered  ac  each  of  the  vortices.  The 
circular  streamlines  are  the  surfaces  of  the  cylinders  which, 
in  this  two-dimensional  framework,  represent  the  bubbles. 


2.2  A  Simple  Hollow  Vortex  Sheet 

Sir  Horace  Lamb  [19]  reports  that  in  an  infinite  expanse  of 

fluid  a  vortex  sheet  moves  with  a  velocity  equal  to  the  average 

of  the  velocities  above  and  below  the  vortex  sheet  U  =  [u^  +  M^ll. 

Furthermore,  in  [13]  it  is  shown  that  a  plane  infinite  row  of 

vortices  does  not  induce  a  velocity  in  itself.  Similarly,  one  can 

show  that  a  plane  infinite  row  of  dipoles  induces  a  velocity  field 

2  2 

in  itself  given  by  Dtt  /3X  ,  where  D  is  the  dipole  strength  and  X  is 

Che  spacing  between  the  dipoles. 

The  circulation  of  the  vortices  is  uniquely  determined  by  the 

average  velocity  across  the  vortex  sheet  U  as  well  as  the  spacing 

of  the  vortices,  X.  From  symmetry  considerations,  ail  Che  dipoles 

in  the  undisturbed  state  must  have  the  same  strength  and 

orientation.  Since  the  vortices  are  being  displaced  with  a 

velocity  U,  and  Che  dipoles  create  a  velocity  field  given  by 
2  3 

Dt7  /3X  ,  Che  normal  velocity  condition  at  each  hollow  vortex  is 
satisfied  if. 


(2.2.1) 


where  R  =  a/X. 
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We  notice  that  the  dipole  strenL.ch  will  be  zero  if  either  R 
or  U  are  zero.  In  the  following  analysis,  we  assume  that  the 
small  disturbance  does  not  affect  the  dipole  strength  and 
orientation  in  a  significant  way.  This  assumption  is  acceptable 
because  we  are  only  Interested  in  small  disturbances. 

According  to  Figure  2.1,  the  original  position  of  the  hollow 
vortices  is  given  by  (mX,0)  while  the  disturbed  state  is 
identified  by  fx  +  mX,  ym),  where  m  is  an  integer.  We  now  focus 
our  attention  on  a  typical  vortex,  say  investigate  how 

the  velocity  field  introduced  by  the  disturbance  affects  the 
stability. 

It  can  be  easily  shown  that  the  velocity  field  created  at  a 
particular  point,  say  (^g>yQ)  displaced  system  of  dipoles 

and  vortices  is  given  by 


(2.2.2a) 


X  -  X 

-  Y  Y  _o _ m  _  r  j_ 

,  ,2  i  2  o  ,2  ^  m 


(2.2.2b) 


(2.2.2c) 


and 
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2D  .  \ 


(2.2.3b) 


where,  in  these  expressions,  only  the  first  order  terms  in  the 
displacements  have  been  kept.  Putting  these  contributions  together, 
one  finds  the  velocity  components  at  the  hollow  vortex  l^g»yQ) 


y  -  y  ^  X  -  X 
-  Y  r  _o _ m  2^  r  o  m 

2  A  2  ^  3  ’ 

27rA  „  m  X  m  m 


(2.2.4a) 


dt  0  ,2  ^  2 

^  TT  A  ^  rn 


.  ^o  '  \  2D  ^o  ~ 
L  0  O  2,  -5 


,3  ^  3 

^  m 


(2.2.4b) 


where  our  reference  system  moves  with  the  rigid  body  velocity  of  the 
vortices  and  dipoles. 

Following  von  Karman,  we  consider  disturbances  of  the  type. 


ime  ime 

x=xe  :  y=ye 
mo  mo 


(2.2.5) 


where  0  <  0  <  27r.  With  this  the  system  of  equations  can  be 
transformed  into 


dx  yy  2Dx 

-  ^A^(e)  --f  iA,(9) 
2ttX  X 


(2.2.6a) 


rT".' j-.  j-  -■ 


d  y  iX  2Dy 

^  - - ^  A  (  9)  +  — iA,  ( 0 ) 


(2.2.6b) 


where  we  have  made  use  of  the  following  series  function 


representations : 


9,0  ON  rl-  cos  me 

A  (6)  =  ^  (2^  -  0)  =  )  - - - 

V  2  '“2 

m 


(2.2.7a) 


Aj(e) 


o  o  ,  o  2^  V  sin  m9 
3710  +  2ir  J  =  }_  - r — 


(2.2.7b) 


Transcribing  the  system  of  equations  to  matrix  notation,  i.e., 


[X]  =  [C][X] 


(2.2.8) 


and  assuming  its  solution  to  be  of  the  form, 


[X]  =  [Y]e' 


(2.2.9) 


one  finds  that  the  problem  is  reduced  to  finding  the  eigenvalues  of 


the  matrix  [C].  These  are 


-)  -  (-r) 


(2.2.10) 


As  expected,  we  see  that  in  the  limit  of  vanishing  core  radius 


we  recover  von  Karraan's  result  of  unconditional  instabilitv  of  a 


simple  vortex  sheet.  If  we  define  a  nondimens ional  parameter  given 


by  C  =  Y1/  277D,  then  the  simple  hollow  vortex  sheet  will  be  stable  if 
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<  0  ,  (2.2.11) 

V 

is  satisfied  for  all  frequencies  of  the  disturbance.  We  see  that 
for  9  =  TT  the  stability  interval  is  reduced  to  zero,  thus  there 
is  no  range  of  combinations  of  circulation,  vortex  spacing  or 
dipole  strength  that  can  make  the  hollow  vortex  sheet  stable. 

2 . 3  A  Symmetric  Hollow  Vortex  Sheet 

In  the  realm  of  potential  flow  theory,  another  possible 
explanation  for  the  apparent  stability  of  the  bubble  layer  is 
presence  of  the  wall.  With  the  same  underlying  assumptions  as 
in  the  last  section,  we  attempt  to  consider  the  presence  of  the 
wall.  The  wall  is  located  at  a  distance  b/2  from  the  initially 
undisturbed  hollow  vortex  sheet.  The  zero  normal  wall  velocity 
boundary  condition  is  satisfied  by  considering  the  image  system 
as  shown  in  Figure  2.2. 

In  the  same  way  that  a  vortex  sheet  in  the  presence  of  a  wall 

induces  a  velocity  in  itself  [13],  so  does  a  dipole  sheet. 

Consequently,  the  hollow  vortices  in  the  disturbed  system  will 

be  identified  by  fmX  +  (U  +  U,)t  +  x  ,  b/2  +  v  '  for  the  upper 

V  d  m  '  m  ' 

row  and  fnA  +  (U  +  U, )t  +  x  ,  -  b/2  +  v  for  the  lower  row, 

^  V  d  m  '  m'  ’ 

where  U  and  U,  are  the  velocities  induced  bv  the  vortex  and 
V  d 

dipole  systems  in  themselves,  respectively. 


ooooo 


Directing  our  attention  to  the  hollow  vortex  ;x  ,v  I  in  the 

'  o  ■  o 


upper  row  of  hollow  vortices,  one  can  show  after  some  algebra, 
that  the  components  of  the  velocity  field  created  by  the  lower 
row  of  vortices  and  dipoles  are  given  by 


r  2  2  2-, 

y  b  Yr  nX-b  . 

1  - -  _  _ —  I  1 - L—  ,  V  -  V  t 

TO  n^ 

"  (n"X“  +  b“]  ^  {nX''  +  b") 


-  Y  '■  2n  Xb  (  \ 

—  ‘  — - -2  (’■o  -  "n)  ■ 

"  (n^X^  +  b“) 


(2.3. la) 


2  2  2 

Y  -  b^  , 


dt  2it  - 

'v  n  2  2,2  .2', 

"  [n  X  +  b  J 


tr  IX  -  X 

2^0  n  ■ 


^  h  I 


-  ^  - - ^  ^y^  -  y^_ 

n  (n^-  +  b"] 


(2.3.1b) 


dx  ] 
0  I 

.  —  n  V 

0  0  O 

n^X  - 

,,  r  2nX(n'"x^ 

dt  ' 

~  ^  L 

a  ^  ^  i 

^  n 

o  0 

f  n  “  X  +  b  "  ? 

"  (n^X^  + 

X  -  X  j 

3  o  n ' 


^  b  ,  30“  X-  -  b‘-  I 
-  2D  ^  ,v  -  V  I 

2  n  o  3  ^  o 

[n  X"  +  b“ j 


(:.3.:a' 


n 


and 
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d  V 


dt 


.nv,2  2.2^ 

-_  r  n Al  3b  -  n  A  !  / 

=  2D  )  - - - 7--  (y  -Vi 

^  OT  o^'-o  -n' 

u  n  r . ^ 


"  (n  A  +  b"] 


-  2D  [ 


.  2  2 
bf3n  A 


(n^A^  +  b^] 


(x  -  X  1 

^  A  n  ^ 


(2.3.2b) 


respectively . 

In  looking  at  our  typical  vortex,  which  is  located  in  the 
upper  row,  we  also  have  to  take  into  account  the  contribution  of 
the  upper  row  to  the  flowfield.  This  contribution  is  expressed 
by  Eq.  (2.2.2)  and  Eq .  (2.2.3),  where  the  vortices  have  the 
opposite  orientation  in  order  to  satisfy  the  zero  normal  wall 
velocity  boundary  condition.  Thus,  the  equations  for  the  velocity 
components  of  the  hollow  vortex  identified  by  become 


dx 

c 

dt 


y  -  y  X  -  X  dx 

Y  '-  o  m  2DrO  m  o 

-I-y  ^  - . -  +  V  -  *  - 


2  TT  X  IT]  ^ 


^  m  "> 


dt 


dx 

o  I 

d  '  V 


,  (2.3.3a) 


and 


dy  x  -  X  or,  y  “  y  dy 

_ o_  _  Y  y  _o _ m  2D  Y  o _ m  o 

2  ~  3  2  ' 


dt 


.2  ^ 
m 


m 
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dt 


dy^ 

dt 


,  (2.3.3b) 


where  the  rigid  body  contributions  from  the  vortices  and  dipoles 
have  been  subtracted.  The  rigid  body  velocity  of  the  system  is 
given  by 
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2  2 
Dtt  X  ,  r-rrbx  ti  D  ,  2  ,-ub-, 

77 "  2i  - -r '  ’ 

j  A  A 


where,  as  can  be  easily  shown  [14], 


coth  f^l 


2tt  2,2  ,.2 

n  n  X  +  b 


2  2  2  2 
■trD  ,2/Trb',  nX  -b 

—  csch  [— )  =  D  I  - : 

^  "  [n^X^  ,  7)' 


Again,  following  von  Karman,  we  define  the  following 


disturbances , 
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x=xii  ;y=yJi 
mo  m  o 


-  „in9  -  ind 

X  =  X  I  ;  v=yll 

n  o  'no 


■-  . 


(2.3.4) 


(2.3.5a) 


(2.3.5b) 


(2.3.6) 


for  the  upper  and  lower  rows,  respectively. 

With  the  function  representations  of  the  resulting  series 
(see  Appendix  A),  we  find 


=  y„  {-7r  [a  (6)  -  c  (0)1  -  ^  C.(0)}  -  LX  ^  A,(e) 
dt  o.,2‘'V  V  •'  ,3d  0,3d 

2itX  X  X 


-*•  y  c  (9)  +  —  C.(9)1 

°  2.x-  "  X^  ' 


B  (9)  +  —  B.(0)  ;  , 

V  ,3d 


(2.3.7a) 


VV 


C  (6)  =  - 


TT-coshka  Tt6sinhk(Tt  -  ■3) 


s  inh  “ktr 


sinhku 


.2  ,2 


(2.3.8d) 


and  k  =  b/X. 


As  in  the  problem  of  a  simple  symmetric  vortex  sheet  there  are 
two  types  of  solutions  for  the  system  of  equations  (2.2.7),  i.e., 


[*29>  -  <:„«»  -  C„<9)]  -  ^  *  Ch»)]) 


9'  °  2.x9  ” 


+  X  i  B  (6)  +  —  [B.(e)  -  2A.(e)]}  , 

°  2trx2  ‘I  ^ 


(2.3.9a) 


=  K(9>  -  c„(0)  +  cje)]  -  ^  [C^(0)  -  C^(6)]} 


°  2trx2 


+  y^i  B^(e)  +  ^  [2A^(8)  +  B^(0)]}  ,  (2.3.9b) 


2^x2 


for  X  =  X  and  y  =  -  v  and  , 
0  0  o  o 


77^  =  y„  K(9)  -  c,(n)  +  c„(6)]  -  -  [c^(n)  -  c,(e) 


2nx2  ^ 


.  i  iB.O)  +  2A^(0)]  +  ^-7  .  (2.3.10a) 


-  .2  V 
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dt 


X 

o 


ZirX' 


[a  (0) 


C  (0)  -  C  (e)j  -  ^  'c,(0)  +  C,(9)'} 

\T  \r  -i  S  •  H  n  - 


+  y  i  i2A,(6)  -  B.(0)1  -  B  f9)}  ,  (2.3.10b) 

o  ,  3  d  d  -  2  V 

X  2ttX 

for  X  =  -  X  and  v  =  y  . 

0  o  '00 

We  realize  immediately  that  if  we  take  away  the  dipole 
contributions  the  two  systems  above  reduce  to  those  of 
von  Karman's  as  reported  in  [13].  Each  of  the  two  systems  of 
equations  can  be  written  in  matrix  form  as  represented  by 
Eq.  (2.2.8).  If  a  solution  of  the  form  of  Eq .  (2.2.9)  is  also 
assumed,  the  stability  condition,  which  turns  out  to  be  the 
same  for  the  systems  (2.3.9)  and  (2.3.10),  can  be  easily  found 
to  be 


{[a^(9)  -  C^(0)]^  -  C^(e)n  -  2c  {C^(9)C^(9) 

+  C^(0)  [A^(9)  -  C^(0)]}  +  [C^(0)^  -  C^(0)^  -  ^A^(9)“’  «  0  . 

(2.3.11) 

where  c  =  yX/2ttD. 

The  stability  of  the  system,  which  has  to  be  satisfied  for 
all  0e[O,2TT)  depends  on  the  initial  wall  distance  of  the  vortices 
(b/2)  as  well  as  the  wavelength  of  the  disturbance  (X).  This 
dependency,  expressed  by  the  dimensionless  parameter  k,  determines 
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the  stability  envelope  given  in  Figure  2.3.  The  shape  or  the 
Stability  curves  sugges^rs  that  the  presence  of  the  wall  has  a  ver 
strong  stabilizing  influence  on  the  hollow  vortex  sheet.  '•■■'e  also 
note  that  the  result  of  unconditional  instabilitv  for  a  hollow 
vortex  sheet,  obtained  in  the  last  section,  is  recovered  just  a 
few  wavelengths  away  from  the  wall. 
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CHAPTER  III 

EVOLUTION  OF  A  HOLLOW  VORTEX  SHEET 

3. 1  Introduction 

The  results  obtained  in  the  last  section  can  be  observed 
qualitatively  from  the  numerical  evolution  of  the  hollow  vortex 
sheet  submitted  to  an  initial  disturbance.  The  numerical  results 

■m 

stem  from  the  consideration  of  the  full  nonlinear  eauations  of 
motion.  Even  though  more  complex,  the  following  development  has 
the  advantage  of  looking  at  the  problem  from  a  dynamical  point  of 
view. 

The  solution  presented  below  follows  the  conceptual  approach 
and  the  steps  of  D.  Rosenhead  [12),  who  by  making  use  of 
von  Karman's  analysis,  tried  to  compute  for  the  first  time  the 
evolution  of  a  vortex  sheet.  It  is  known  today  [10]  that 
Rosenhead 's  computations  were  incorrect.  Indeed,  some  authors 
believe  that  the  discrete  vortex  representation  of  a  surface 
velocity  discontinuity  may  not  be  appropriate.  Nevertheless,  this 
has  been  the  approach  used  by  C.  Y.  Chow  [11]  and  R.  F.  Hama, 
et  al.  [10]  to  compute  the  evolution  of  the  vortex  sheet.  The 
latter  authors  point  out  interesting  first  order  vorticity  effects 
that  had  been  disregarded  in  all  previous  analyses  and  which 
cannot  be  accounted  for  in  the  vortex  representation  of  the 
surface  velocity  discontinuity.  Through  a  simple  model  they  also 
looked  at  how  the  background  vorticity  affects  the  evolution  of 


the  vortex  sheet. 
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The  analysis  presented  below  is  suggested  in  the  first  chapter 
of  a  most  interesting  dissertation  by  G.  R.  Baker  [15]  on  vortex 
dynamics . 

As  in  the  last  chapter,  we  attempt  to  model  the  hollow  vortices 
by  using  dipoles.  In  contrast  with  Chapter  II,  however,  the 
orientation  and  strength  of  the  dipoles  depends  on  the  character¬ 
istics  of  the  flowfield  past  each  of  the  hollow  vortices,  vide 
Figure  3.1. 

There  are  two  characteristic  lengths  in  the  problem,  i.e.,  the 
cylinder  radius  (a)  and  the  wavelength  (X).  Since  the  ratio  of  the 
hollow  core  radius  to  the  wavelength  (a  =  ira/X)  is  very  small,  we 
will  take  into  account  only  the  first  order  contribution  of  the 
disturbance  to  the  solution.  From  Helmholtz's  theorems  it  follows 
that  to  solve  the  point  vortex  problem  one  needs  only  to  take  into 
account  the  kinematics  of  the  flowfield.  To  solve  the  problem  of 
the  hollow  vortices,  however,  one  must  consider  both  the  kinematics 
and  the  dynamics  of  the  system  of  hollow  vortices.  The  forces 
created  by  the  flowfield  on  each  of  the  hollow  vortices  will 
determine  their  evolution  which  in  turn  determines  the  dipole 
strength  and  orientation. 

3.2  The  Dipole  Strength  of  the  Svstem 

Since  the  flowfield  under  consideration  is  a  two-dimensional 
inviscid  and  irrotational  flow,  one  can  define  both  a  stream- 
function  (fi)  and  a  velocity  potential  (li).  The  complex  potential 
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and  the  complex  velocity  Is  obtained  by  taking  the  derivative  of 
W(z)  with  respect  to  z.  The  velocity  components  then  relate  to 
the  streamf unction  and  the  velocity  potential  by  the  Cauchy- 
Riemann  conditions. 

It  can  be  easily  shown  [16]  that  the  complex  potentials  of 
a  dipole  and  a  vortex,  located  at  z^,  are  given  by 


W(z)=D/fz-zl  (3.3.2a) 

Q  O 

and 

W  (z)  =  ir  logfz  -  z  ]  (3.3.2b) 

respectively.  In  Eq .  (3.3.2b),  F  =  Y/2Tr  and  the  positive  sense  is 

taken  as  the  counterclockwise.  From  expression  (3.3.2a)  we  notice 

2  2 

that  the  dipole  strength  D  is  proportional  to  a  ,  i.e.,  D  =  La  . 

The  complex  potential  of  a  plane  infinite  row  of  vortices 
can  be  shown  [16]  to  be  given  analytically  by 

(z)  =  ir  log  [sin  y  (z  -  z^)]  (3.3.3) 


where  z^  represents  the  position  of  the  vortex  closest  to  the 
origin  and  A  the  distance  between  vortices,  vide  Figure  2.1. 

From  the  application  of  the  Circulation  Theorem  over  one  wave¬ 
length  of  the  vortex  sheet,  we  find  that  F  =  UX/ttN,  where  N  is 
the  number  of  hollow  vortices  contained  in  one  wavelength  of  the 
disturbance . 

If  we  differentiate  Eq .  (2.3.2b)  and  compare  it  with 
Eq .  (2.3.2a),  we  conclude  that 


(3.3.4) 
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Since  the  comp]  ,'.x  velocity  potential  of  an  Infinite  row  of 
vortices  is  obtained  by  superposing  the  velocity  potentials  of 
each  of  the  vortices,  we  can  write 


«d  = 


iP 

r 


dW  (z) 

V 


DtT  ,  Tt  ^  N 

-cot  -  (z  -  zj 


dz 


(3.3.5) 


Thus,  the  complex  velocity  fields  created  by  a  plane  infinite 
row  of  dipoles  and  vortices  are  given  in  closed  analytic  form  by 


dW 


«  i  Ftt 


1  1  Tf  _  ^  r-  'I 

=  nr  T 


(3.3.6a) 


dW, 


dz 


Dir  2  TT  / 

—  =5=  I  (»  -  . 


(3.3.6b) 


respectively. 

Since  we  are  working  with  a  potential  flowfield,  we  are 
ultimately  solving  Laplace's  equation,  for  which  both  the  real 
and  imaginary  parts  of  each  of  the  above  complex  potentials  are 
solutions.  In  this  case,  the  superposition  principle  is  valid 
and  the  complex  potential  of  the  whole  flowfield  will  be  composed 
of  contributions  of  the  infinite  rows  of  dipoles  and  vortices  as 
well  as  their  images  in  each  of  the  hollow  vortices. 

As  in  the  case  of  the  simple  point  vortex  system,  the 
resulting  equations  are  nonlinear.  In  order  to  keep  the 
analysis  within  bounds,  we  are  forced  to  consider  only  the 
leading  contribution  of  the  hollow  vortices  to  the  flowfield. 

To  evaluate  the  order  of  magnitude  of  the  contribution  of  the 
image  system,  we  invoke  Milne  and  Thomson's  Circle  Theorem 
[17]  and  look  at  the  image  system  of  a  dipole  and  a  vortex  in 
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the  presence  of  a  cylinder*  In  the  two  simple  cises  of  a 


cylinder  located  at  the  origin  and  a  dipole  or  a  vortex  located 


at  2  ,  we  find  that  the  complex  velocities  can  be  expressed  as 
o 


d  [z-zi  zz  zz 

^  o  •'  o  o 


(3.3.7a) 


dz  I  z  -  2 

'  V  ^  o  ^  z  z 
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r  1  +  ^  + 

')  L  ^ 


(l.i.lh) 


respectively. 


Since  Che  dipole  strengths  (D)  are  of  order  a',  we  conclude 


Chat  to  this  order  of  approximation  the  image  system  due  to  the 


dipoles  can  be  neglected  while  that  of  the  vortices  must  be 


considered.  Also,  to  this  order  of  magnitude,  we  note  that  to  cake 


into  account  the  contribution  of  the  vortex  image  system  is  the 


same  as  to  consider  a  dipole  at  the  origin  with  a  strength  given 


by  ifa^/z 


‘^ince  Che  flow  under  consideration  satisfies  the  premises  of 


the  Circle  Theorem,  the  complex  potential  of  the  image  system  of 


a  plane  row  of  vortices  in  a  cylinder  located  at  z  is  given  bv 

o' 


W  (z)  =  ir  —7^  F  [z  !  cot  —  (z  -  z  1 
V,  ,2  k  ^  o '  X  ^  o' 

k  X 


(3.3.«) 


where  the  function  F  [z  ]  depends  on  the  position  of  all  the  ocher 


vortices  and  Che  index  k  denotes  the  row  to  which  these  vortices 


belong.  It  can  be  shown  that  [lA], 


=  I  I - ^ 


=  IT  cot  Y  (z  -  zj 


(3.3.9) 
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With  the  above  in  mind,  we  recognize  f jur  major  contributions 
to  the  flowfield  complex  potential: 

•  The  complex  potential  of  the  infinite  row  of  vortices. 

•  The  complex  potential  of  the  image  system  of  the 

vortices  in  each  of  the  hollow  vortices. 

®  The  complex  potential  of  the  infinite  row  of  dipoles. 

•  The  complex  potential  of  the  flowfield  due  to  the 

motion  of  the  hollow  vortices. 

The  last  contribution  is  equivalent  to  the  action  in  each  of 
the  hollow  vortices  of  another  dipole  the  strength  of  which  depends 
on  the  velocity  and  the  diameter  of  the  hollow  vortex  and  is  given 


D  =  Ua 


(3.3. 10) 


where  U  is  a  complex  number  representing  the  velocity  of  the  hollc 


vortex. 


In  the  realm  of  potential  flows,  the  complex  potential  of  the 
whole  flow  field  becomes 


W(z)  =  I  {ir  log[sin  y  (z  “  ^ 


+  G  cot  Y  (z  -  z  )} 

n  X  X  ^  n-*^ 


(3.3.11) 


where , 


L  +  U 
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(3.3.12) 
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where  Che  last  term  represents  the  velocity  field  created  at  the 
hollow  vortex  by  the  dipoles  in  its  own  row. 

In  order  to  satisfy  the  zero  normal  velocity  boundary 
condition  at  Che  surface  of  each  of  the  hollow  vortices,  we 
must  require  that 

-  L  =  0  m  =  .  (3.3.14) 

m 
m 


3z 


This  expression  represents  a  linear  system  of  N  equations 
and  N  unknowns,  which  enables  the  computation,  at  each  time  step, 
of  the  dipole  strengths  of  Che  system  of  hollow  vortices.  For  a 
better  understanding  vide  Figure  3.1  where  the  various  quantities 
involved  in  the  computation  of  the  dipole  strengths  are  shown 
schematically.  To  a  first  order  of  magnitude,  the  above  equation 
reduces  to 


N 

L  =  cot  "7  (2  “  z  ]  ®  =  1,...,N  .  (3.3.15) 

m  X  ^  A  ^  m  n-* 

n  =  l 

n^ 

Thus,  to  a  first  order  of  magnitude,  the  dipole  strengths  are 
determined  by  the  vortex  flowfield.  This  is  to  be  expected  since 


in  a  two-dimensional  flowfield  the  vortex  velocity  field  vanishes 
as  1/z  while  the  dipole  velocity  field  vanishes  as  l/z*". 


(3.3. 16) 
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The  rate  of  chanj^e  of  the  dipole  strengths  is  given  by 


dL  .  ^  2  ^ 

m  liTT  ^  r  \  2  "n  r 

T—  - - (u  -  U  )  CSC  -r  ~  ^  ■  > 

at  2  n'  n' 

^  n=  1 


where  again  m  =  1,...,N. 


3. 3  The  Dynamics  of  the  System 

In  order  to  track  the  motion  of  the  hollow  vortices,  we  look 
at  their  equations  of  motion.  In  complex  variables,  this 
corresponds  to  the  application  of  Blasius  Theorem  [16].  It  can 
be  easily  shown  that,  if  one  neglects  the  rotation  of  the 
boundary,  the  generalized  form  of  Blasius  Theorem  reduces  to 


X  -  iY 


.  r  ^  1  .  f 

■  *  2  /  fe)  '1' 

ra  m 


^  ‘“l  /  sT *  V  SF-  ,(3.3.17) 

ra 

where  the  left  hand  side  is  the  inertia  term  in  aopropriate  complex 
notation  and  the  right  hand  side  arises  from  the  integration  of 
Euler's  equations  around  a  closed  contour  enclosing  the  body's 
cross  section.  In  Eq .  (3.3.17)  represents  the  contour  of  the 
cross  section  of  the  m^^  hollow  vortex.  As  is  shown  in  Appendix  B, 
each  of  the  five  terms  of  the  above  equation  is  found  to  give 
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(3.3.  18c) 
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(3.3.18d) 
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X  -  iY|  =  wp  7~ 
'  m  g  d  C 


(3.3. 18e) 


where  ra  =  1 , . . . , N. 

Equations  (3.3.18a)  and  (3.3. 18d)  account  for  the  unsteadiness 
of  the  flow  and  are  thus  equivalent  to  the  added  or  virtual  mass 
term.  Equation  (3.3.18b)  gives  both  the  first  and  second  order 
steady  state  flow  contributions.  In  Eq.  (3.3. I8c),  we  recognize 
Kutta-Joukowski 's  Theorem  [15]. 

In  the  limit  of  vanishing  core  radius,  we  find 
N 

U  =-irT  y  cot^fz  -zl  m=l,...,N  ,  (3.3.19) 
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which  are  exactly  the  equations  of  motion  of  the  point  vortices 
in  a  sheet  as  given  by  Sir  Horace  Lamb  [13]. 

In  dimensionless  variables,  the  equations  of  motion  reduce  to 
a  system  of  first  order  ordinary  differential  equations  of  the 


»  dU 

a  =  S  (U,Z) 

d  t  m - 


m -  =  U 

dt  m 


m  !,•••, N  , 


(3.3.20) 


where  IJ  and  ^  represent  the  velocity  and  position  arrays  of  the 
system  of  hollow  vortices.  In  arriving  at  the  system  of 
Eqs.  (3.3.20),  we  used  the  fact  that  Pg/pi  <<  1.  After  some 
algebra,  the  source  term  Sju  is  found  to  be  given  by 


S  (U,Z)  =  -  {i  +  2Im}f  T  [U  -  U  Icsc^irfz  -  z  )1 

m -  N  ‘■^^m  n'  m^'* 
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So,  from  Eq.  (3.3.21),  we  recognize  tbac  in  Che  syscein  (3.3.20) 
there  is  actually  imbedded  the  zeroth  order  problem  that  represents 
Che  motion  of  the  point  vortices. 

In  order  to  solve  the  system  (3.3.20),  given  the  initial 
position  of  the  hollow  vortices  the  zeroth  order  problem, 
represented  by  Eq .  (3.3.19),  determines  the  initial  velocity  of 
the  hollow  vortices.  With  the  initial  vectors  IJ  and  the 
corrections  to  the  velocities  of  the  system  of  vortices  are 
computed  from 
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dt  m - 
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(3.3.23) 

Integrating  the  second  equation  of  the  system  (3.3.20)  and 
using  as  the  source  term  the  result  obtained  from  (3.3.22)  one 
obtains  the  corrections  for  the  the  positions  of  system  of  hollow 
vortices.  Thus,  the  new  position  for  the  hollow  vortex  sheet  is 
described  by 


z  =  z  +  a  z 
m  oni  cm 


in  *  y 


(3.3.24) 


where  is  obtained  as  described  above  and  Zq^i  is  the  result  of 

the  integration  of  (3.3.19)  and  illustrated  for  one  case  in 
Figure  3.2. 
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A  fourth  order  Runge-Kutca  solver  was  used  to  perform  the 
computations  described  above.  A  flow  chart  as  well  as  the  code 
are  included  in  Appendix  C.  Supporting  our  previous  findings, 
comparison  of  Figures  3.2,  3.4  and  3.5  suggests  that  the 
consideration  of  the  velocity  field  due  to  the  dipoles  is  a 
destabilizing  factor.  The  results  indicate  also  that  the  larger 
the  cylinder  radius,  the  shorter  the  roll  up  time.  In  the  two 
cases  presented,  the  cylinder  diameters  make  up  25%  and  46%  of 
the  wavelength. 

3.4  The  Image  System 

Since  the  small  disturbance  analysis  carried  out  in  Chapter  II 
suggests  that  the  presence  of  the  wall  is  the  stabilizing  factor 
of  the  hollow  vortex  sheet,  we  introduce  an  image  system  as 
illustrated  in  Figure  2.2. 

With  the  complex  potentials  of  a  dipole  and  a  vortex,  located 
at  Zq,  given  by  Eqs.  (3.3.2),  it  is  clear  that  the  corresponding 
complex  potential  of  the  image  system  in  a  horizontal  plane  is 
given  by 

W(z)=D/fz-z]  ,  (3.4.1a) 

a  ^  ^  o 

and 

W  (z)  =  -  ir  log(z  -  z  )  ,  (3.4.1b) 

Vi  o 

respectively.  The  bars,  of  course,  mean  the  complex  conjugate 
of  the  quantity. 

With  an  analysis  similar  to  that  of  the  previous  sections,  one 
concludes  that  the  complex  potential  of  the  hollow  vortex  sheet 
with  the  corresponding  image  system  is  given  by 


W(z) 


=  \  {ir  log[sin 

n  =  l 


9 


X  (==  -  -  ir  log[sin  X  (z  - 

■  "n)  i 


where 


G' 

n 


N 

L'  +  U  -  i  I  T 

k=l 

k^n 


-  z. 


.  r  TT  ^  IT  r- 

+  1  —  )  cot  —  [z  -  Z, 

X  ^  X  '■  n 

k  =  l 
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In  the  last  expression,  the  last  term  stems  from  the  image  in 
each  of  the  vortices  of  the  wall  image  system. 

Enforcing  the  zero  normal  velocity  boundary  condition  at  the 
surface  of  each  of  the  hollow  vortices,  as  expressed  by 
Eq .  (3.3.14),  we  find  that  for  the  symmetric  hollow  vortex  sheet, 
to  the  first  order  magnitude,  we  have 
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(3.4.4) 


With  an  argument  identical  to  that  of  the  last  section,  the 


motion  of  the  hollow  vortices  in  the  presence  of  the  wall  is  still 
governed  by  Eq .  (3.3.17).  Upon  integration  with  the  complex 


potenti";!  as  defined  by  Eq .  (3.4.2),  we  conclude  that  only  the  term 
corresponding  to  Eq .  (3.3.18b)  is  modified  in  the  equations  of 
motion,  i.e.. 
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The  other  two  integrals  are  not  affected  because  the  new 
functional  terms  added  to  the  complex  potential  are  always 
analytic  Inside  the  contour  of  integration  and  consequently 
vanish.  While  in  the  case  of  the  hollow  vortex  sheet  in  an 
infinite  expanse  of  fluid,  the  velocity  scale  is  the  farstream 
induced  velocity  (U),  in  this  case  it  is  the  rigid  body  velocity 
of  the  hollow  vortex  sheet  (V).  Thus,  the  solution  of  the  problem 
with  the  image  system  is  still  obtained  through  the  system 
(3.3.20),  where  now  the  source  term  is  given  by 
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,N.  In  the  above  equation,  u  =  nr/XV. 


The  zeroth  order  of  magnitude  problem  is,  of  course,  dictated 
by  the  complex  potential  of  the  symmetric  arrangement  of  vortices. 
This  complex  potential  falls  out  of  Eq .  (3.4.2)  in  the  limit  of 
vanishing  hollow  vortex  radius.  The  corresponding  complex  velocity 
field  is  given  by 


U  =  -  ioj  {  y  cot  irfz  -  z  )  -  y  cot  mfz  -  z  1} 
m  ^mn'^  *-  -mn^ 

n=l  n=l 


m  *  I,...,N  , 


(3.4.7) 


where  the  indices  refer  to  the  vortices  in  the  upper  row. 

The  integration  procedure  for  this  problem  is  entirely 
identical  to  that  described  in  the  previous  section.  For 
computational  convenience,  but  certainly  to  the  detriment  of  the 
CPU  time  required,  we  used  in  both  instances  complex  notation 


in  the  course  of  the  computations. 

The  stability  range  of  the  symmetric  hollow  vortex  sheet  found 
in  the  last  chapter  imposes  restrictions  chat  cannot  be  satisfied  if 


the  physical  dimensions  of  the  cylinders  are  taken  into  account. 

In  order  to  satisfy  these  impositions,  we  have  to  consider  a  very 
large  number  of  vortices,  which  renders  the  computations 
impractical.  According  to  Figure  2.3  a  stable  vortex  hollow  vortex 
sheet  can  be  found  farther  away  from  the  wall  if  the  hollow  vortex 
circulation  is  small  enough.  Once  more  the  computations  are 
impractical  because  the  hollow  cores  only  introduce  a  second  order 
correction  to  the  point  vortex  motion.  As  a  consequence,  we  were 
unable  to  corroborate  numerically  the  results  found  in  Section  2.3. 
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The  small  size  of  the  bubbles  in  the  boundary  layer  suggests 
that  they  are  endowed  with  a  very  small  rotation.  From  the  previous 
analysis  the  small  circulation  along  with  the  presence  of  the  wall 
are  the  key  reasons  for  the  apparent  stability  of  the  bubbles  in  the 
shear  flow.  Other  contributing  factors,  which  this  analysis 
completely  disregards,  are  the  effect  of  viscosity  which  may  be 
stabilizing  and  the  presence  of  other  bubble  layers. 

For  illustration  purposes,  the  evolution  of  a  point  vortex 
sheet  in  the  presence  of  a  wall  is  given  in  Figure  3.3.  The 
pattern,  clearly  unstable,  in  support  of  von  Karman's  conclusion 
[13],  and  is  suggestive  of  the  mechanism  of  drag  reduction  which 
we  discuss  in  the  next  chapter. 
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CHAPTER  IV 


THE  DRAG  REDUCTION 


4. 1  Introduction 

As  mentioned  in  Chapter  I,  the  promise  of  substantial  skin 
friction  reduction  through  boundary  layer  microbubble  injection 
excited  the  interest  of  many  engineers.  The  obvious  potential 
applications  of  this  concept  of  drag  reduction  led  some  of  these 
researchers  to  attempt  to  model  the  boundary  layer  flow  with 
raicrobubbles  [18].  It  should  be  noted  that  these  attempts  were  of 
a  probing  nature.  The  aforementioned  researchers  claim  to  have 
only  assessed  the  percentage  of  the  drag  reduction  that  can  be 
attributed  to  a  mean  flow  density  and  viscosity  variation. 

In  the  following,  a  more  comprehensive  attempt  is  presented. 
Like  previous  works,  the  direction  adopted  has  shortcomings  which 
are  imposed  by  the  extreme  complexity  of  the  flow  with  its  poorly 
understood  structure.  Even  though  the  modeling  concept  is 
basically  the  same,  we  think  that  if  a  phenomenlogical  approach 
is  to  be  adopted,  then  it  should  be  built  upon  the  rough  model 
introduced  in  the  following  sections. 

Some  of  the  major  difficulties  that  have  to  be  dealt  with 
in  this  problem  in  order  to  design  a  reliable  code  are: 

•  the  fact  that  we  have  two  phase  flow, 

•  the  large  range  of  variation  of  the  void  fraction 
of  the  second  phase  in  the  flow, 

•  the  interaction  between  the  second  phase  in  the 
flow  and  the  turbulence  structure  of  the  carrier 


a  the  momentum  transfer  structure  between  the  second 


phase  and  the  mean  flow  in  cross  stream. 

The  solution  of  the  above  problems  represents  a  major 
mathematical,  numerical  modeling  and  experimental  effort.  Part 
of  this  drive,  mostly  on  the  experimental  side,  has  been  under¬ 
way  at  The  Pennsylvania  State  University  for  sometime  [6-8] . 


4.2  The  Constitutive  Relation 

M.  Ishil,  in  a  most  interesting  monograph  [19],  discusses 
from  first  principles  turbulent  two  phase  flows  and  the 
mathematical  tools  within  reach  to  describe  them.  From  the 


complexity  of  the  equations  obtained,  it  is  obvious,  at  least 
for  the  near  future,  that  any  hope  of  modeling  these  flows  lies 
in  the  possibility  of  representing  them  as  some  kind  of 
continuum  single  phase  flow. 

This  being  the  only  path  open,  it  assumes  the  knowledge  of 
the  physical  characteristics  of  the  "fluid".  This  means  that  an 
experimental  effort  must  be  undertaken  in  order  to  study  and 
establish  a  functional  relationship  between  the  viscosity  of  the 
"fluid"  and  the  void  fraction  over  large  ranges.  In  other 
words,  one  must  first  find  the  appropriate  constitutive 
relations  for  the  "fluid". 

To  the  author's  knowledge,  N.  K  Madavan,  et  al .  [18]  in  their 


numerical  study  addressed  this  problem  by  representing  the 
viscosity  of  the  "fluid"  by  Einstein's  relation  [20,21],  i.e. 


u  =  u„(l  +  2.5X) 


(4.2.1) 


as  well  as  Sibree's  viscosity  model  in  order  to  provide  an  upper 
bound  to  the  effects  of  density  and  viscosity  variation.  The 
validity  of  Eq .  (4.2.1),  which  is  obtained  from  theoretical  hydro- 
dynamic  arguments,  is  restricted  to  small  values  of  the  void 
fraction,  X.  Even  in  this  range,  nevertheless,  the  experimental 
work  of  S.  Hinata,  et  al.  [22]  suggests  that  Eq .  (4.2.1)  is  not  a 
good  description  of  the  functional  dependence  of  the  viscosity  on 
the  void  fraction.  The  American  researchers  did  not  seem  to  be 
aware  of  this  work  by  their  Japanese  colleagues.  Indeed  Figure  4.1 
shows  that  Einstein's  model  is  itself  an  upper  bound  fot  the 
viscosity  when  considered  as  a  function  of  the  void  fraction.  For 
this  reason  and  the  fact  that  in  [4]  the  void  fraction  profiles 
are  prescribed  arbitrarily,  the  author  thinks  that  the  conclusions 
in  [4]  should  be  considered  at  most  qualitatively. 

Without  delving  into  the  merits  of  the  experimental  procedure 
used  by  S.  Hinata,  et  al.,  we  note  two  important  facts: 

The  equation  proposed  by  them  as  a  constitutive  relation 
dependent  on  the  void  fraction,  i.e., 

U  =  X(0.45  +  1.3X)Ta~^^^  ,  (4.2.2) 

sp 

where  Ta  =  u^d(du/dy )/a,  is  based  on  void  fractions  that  range  up 
to  0.25.  Even  though  based  on  this  limited  range  of  void  fractions, 
there  is  reason  to  believe  that  Eq .  (4.2.2)  is  valid  in  our  flow. 

Our  argument  is  connected  to  the  coalescence  and  def ormability  of 
the  bubbles.  According  to  [18]  this  can  be  correlated  to  the  non- 
dimensional  number  Ta  with  the  conclusion  that  if  Ta  >  0.15, 
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coalescence  and  def ormabi} Lty  are  prevalent,  Eq.  (4.2.2)  is  not 

valid.  For  a  typical  bubble  in  the  present  flow  [8],  we  find  that 

-2 

the  parameter  Ta  is  of  order  10  •  Some  of  their  experimental 

findings,  as  well  as  the  empirical  relation  (4,2.2),  are  reproduced 
in  Figure  4.1. 

At  lower  void  fractions,  the  empirical  relation  (4.2.2) 
compares  very  favorably  with  Taylor's  expression  [23],  see 
Figure  4.1, 

2 

+  -r  u. 

u  =  2.5X  r — -]  .  (A. 2. 3) 

g  I 

Consequently,  Eq .  (4.2.2)  is  used  in  our  computations,  even 
though  we  have  no  knowledge  of  how  good  it  is  at  the  higher  void 
fractions  of  0.8  reported  by  the  Soviet  researchers  [2]. 

4.3  The  Mechanisms  of  the  Drag  Reduction 

Even  though  the  mechanisms  of  the  microbubble  drag  reduction 
are  still  unclear,  we  present  here  a  chain  of  ideas  that  seems 
to  be  physically  consistent  in  light  of  the  experimental  data 
available.  The  ideas  presented  will  be  used  in  Section  4.5  to 
modify  the  formulation  of  the  Van  Driest  wall  function. 

M.  K.  Madavan,  et  al  [18]  uses,  in  his  numerical  studies, 
as  the  main  mechanism  to  achieve  the  skin  friction  reduction,  a 
modified  Van  Driest  equation  for  the  mixing  length.  The  physical 
motivation  of  his  expression,  however,  is  obscure  and  in  fact 
it  does  not  capture  one  of  the  main  characteristics  of  the  draq 
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reduction  by  microbubble  injection,  i.e.,  the  fact  that  it  is  a 
localized  streamwise  phenomenon.  This  is  a  consequence  of  the 
fact  that  the  void  fraction  is  arbitrarily  prescribed  in  [4] 
and  not  computed. 

The  stability  results  of  the  past  chapters  suggest  and  the 
experimental  results  give  credence  to  the  idea  that  the  draq 
reduction  is  a  result  of  the  dynamic  interaction  of  the  micro¬ 
bubbles  with  the  macroscale  turbulence  structure  of  the  boundarv 
layer. 


Since,  downstream  from  the  place  of  injection,  the  skin 
friction  eventually  recovers  to  its  undisturbed  value,  it  is 
obvious  that  the  mere  presence  of  the  bubbles  in  the  boundary 
layer  is  not  responsible  for  the  drag  reduction.  In  fact,  the 
streamwise  distance  over  which  the  drag  reduction  is  observed 
is  rather  limited  [8].  Furthermore,  the  fact  that  the  skin 
friction  reduction  depends  on  the  volumetric  flow  rate  of  the 
gas  being  injected  supports  the  last  oaragraph  in  two  wavs. 

It  was  observed  that  [6,8): 

o  The  skin  friction  reduction  is  very  small  for  small 
volumetric  flow  rates  of  the  gas  being  injected. 

®  The  skin  friction  reduction  is  partially  lost  if  the 
volumetric  flow  rate  of  the  gas  being  injected  is 


increased  beyond  a  specific  quantity  that  depends 
on  the  free  stream  Reynolds  number.  In  other  words, 
for  a  given  free  stream  Reynolds  number,  there  is 
an  optimum  injection  volume  flow  rate  for  maximum 


skin  friction  reduction. 


I 


i 
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These  two  facts,  seemingly  unrelated,  have  a  common  connection. 
In  the  first  case,  the  lower  bubble  inertia,  the  larger  diameter 
determines  that  the  bubbles  are  forced  to  interchange  momentum  with 
the  smaller  turbulence  eddies  closer  to  the  wall.  This,  allied  to 
the  fact  that  in  this  flow  regime  a  fewer  number  of  bubbles  is 
injected,  determines  that  only  a  small  percentage  of  the  turbulence 
energy  of  the  mean  flow  is  used  to  deflect  the  bubbles  and 
eventually  convect  them.  On  the  other  hand,  the  higher  buoyancy 
of  the  bubbles  tends  to  reduce  the  time  in  thich  they  interact 
with  the  boundary  layer.  The  last  remark  is  supported  in  a  larger 
context  by  N.  K.  Madavan's  experiments  [8]  in  which  two 
gravitational  orientations  of  the  boundary  layer  were  studied.  It 
was  observed  that  when  the  buoyancy  field  acts  against  the  escape  of 
the  bubbles  from  the  boundary  layer  the  bubbles  become  m.ore 
effective  in  reducing  the  skin  friction,  vide  Figures  1.5  and  1.6. 

The  above  ideas  explain  why  an  increase  in  the  injection  flow 
rate  will,  in  general,  lead  to  a  further  decrease  in  the  skin 
friction.  The  higher  inertia  of  the  bubbles  allows  them  to  reach 
larger  eddies  farther  out  in  the  boundary  layer,  while  their 
Smaller  diameter  and  larger  number  makes  them  more  effective 
vehicles  of  momentum  transfer. 

The  fact  that  for  a  given  injection  flow  rate  lower  skin 
friction  reductions  are  observed  for  higher  free  stream  velocities, 
vide  Figure  1.1,  can  also  be  explained  by  the  above  reasoning. 

At  higher  free  stream  velocities,  the  bubbles  are  swept  away  sooner 
by  the  mean  flow  forcing  the  momentum  interchanges  with  the  mean 
flow  to  occur  through  smaller  eddies  closer  to  the  plate.  One 
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should  keep  in  mind  that  the  above  reasoning  assumes  a  zero  pressure 
gradient.  Obviously  the  presence  of  a  pressure  gradient  affects  the 
dynamic  behavior  of  the  bubbles  in  the  boundary  layer.  Consequently 
the  skin  friction  reduction  has  different  characteristics  from  the 
flat  plate  case. 

The  second  experimental  fact  mentioned  above,  however,  seems 
to  contradict  the  preceeding  exposition.  The  answer  to  the 
apparent  contradiction  lies  in  the  structure  of  the  turbulent 
boundary  layer  [24] .  The  turbulent  energy  is  supplied  to  the 
boundary  layer  from  the  free  stream  through  large  scale  structures, 
which  through  vortex  stretching  convey  the  turbulent  energy  to  the 
inner  regions  of  the  boundary  layer. 

With  this  rough  picture  in  mind,  we  conclude  that  for  a  given 
free  stream  flow  regime  the  maximum  skin  friction  reduction  is 
obtained  for  the  injection  volumetric  flow  rate  for  which  the 
bubbles  have  enough  inertia  to  reach  and  interact  with  the  larger 
scale  structures.  In  this  way  the  bubbles  provide  a  breakdown  of 
the  chain  of  turbulent  energy  transmission  to  the  inner  region  of 
the  turbulent  boundary  layer.  The  fact  that  the  skin  friction 
starts  increasing  after  the  injection  flow  rate  has  surpassed  its 
optimum  value  is  clear  because  in  these  circumstances  an  increasing 
number  of  bubbles  has  sufficient  energy  to  bypass  the  large  scale 
structures  without  Interchanging  momentum  with  the  turbulent 
structure  of  the  boundary  layer.  Finally,  the  localized  nature 
of  the  phenomenon  is  also  explained  by  the  fact  that  once  bubbles 
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are  brought  to  the  same  velocity  as  the  mean  flow,  they  possess 
no  more  energy  to  interchange  with  the  flow  and  thus  lose  their 
ability  to  affect  the  flow's  turbulent  structure. 

From  this  picture  one  perceives  four  main  contributors  to  the 
drag  reduction  in  a  boundary  layer  with  microbubbles: 
o  A  decrease  of  the  mean  density  of  the  "fluid", 
o  An  increase  of  the  mean  viscosity  of  the  "fluid", 
o  The  interchange  of  linear  and  angular  momentum 
between  the  bubbles  and  the  turbulence  structure 
of  the  mean  flow,  and 

®  The  added  mass,  or  fluid  entrainment,  of  the 
bubbles  when  they  are  injected. 

At  the  stage  of  current  research  on  this  problem,  it  is 
difficult  to  identify  the  main  contributors,  however,  the  previous 
reasoning  suggests  that  the  last  two  physical  phenomena  are  probablv 
the  most  important. 


4.4  The  System  of  Equations 


The  equations  to  be  integrated  for  laminar  flow  are,  of  course, 
the  steady-state  compressible  Navier-Stok.es  equations  [25],  along 
with  a  transport  equation  for  the  second  phase. 


V  •  (pV)  =  0 

p[v  •  V  v]  =  pB  -  Vp  -  pVx(VxV)  -  7ux(VxVl  , 


V  •  VX  =  0 


as  well  as  their  turbulent  counterparts. 


(4.4.1a) 


(4.4.1b) 


(4.4.1c ) 
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In  Che  system  above,  y  is  the  viscosity  of  Che  mixture  which 
is  related  to  Eo  •  by 


y  =  u  f  1  +  U  } 
sp" 


(4.4.2) 


The  transport  equation  (4.4.1c)  stems  from  the  assumption 
chat  the  bubbles  are  convecced  by  Che  main  flow.  Even  though 
this  is  not  true  at  the  locus  of  iniection,  we  expect  this 
equation  to  capture  some  of  the  kinematics  associated  with  Che 
fact  that  the  bubbles  at  injection  possess  some  velocity  of  their 
own.  Besides  this,  the  adoption  of  Eq .  (4.4.1c)  introduces  in 
Che  system  the  ability  to  compute  the  void  fraction  profile,  as 
opposed  to  its  a  priori  specification  as  is  done  in  [18). 

Equation  (4.4.1c)  can  be  understood  as  a  remanent  from  the 
consideration  of  a  full  two  phase  flow  system.  Clearly  Che 
local  density  in  the  system  (4.4.1)  is  given  by 


p  =  p  (1  -  X)  +  p  X  .  (4.4.3) 

For  a  two-dimensional  flow,  the  system  (4.4.1)  reduces  to 

(4.4.4a) 
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(4.4.4c) 

(4.4.4d) 


where  to  obtain  Eq .  (4.4.4a),  we  used  Fa.  (4.4.4h)  and  Eg.  (4.4.3). 
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Since  we  know  very  little  about  the  momentun  interchanees 
between  the  bubbles  and  the  carrier  phase,  we  will  discard  the 
body  force  field  terms.  In  this  continuum  aoproach,  however,  we 
still  can  take  into  account  the  buoyancy  force  field  of  the  bubbles, 
Neglecting  the  other  dynamic  interactions  may  not  be  a  very  good 
approximation  since  the  experimental  results  summarized  in 
Chapter  I  showed  the  biggest  skin  friction  reduction  close  to  the 
region  where  the  bubbles  are  injected,  a  fact  that  may  be 
associated  with  momentum  interchanges  between  the  bubbles  and  the 
mean  flow.  Dynamic  interactions  between  the  bubbles,  which  are 
probably  important  in  this  region  due  to  the  high  packing  and 
velocity  difference  between  the  two  phases,  had  also  to  be 
disregarded . 

By  virtue  of  the  fact  that 


p^/Pj  «  1  , 


the  system  (4. A. 3)  can  be  further  reduced  to 
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For  this  problem,  the  laminar  flow  regime  develops  in  the 


absence  of  bubbles  and  thus  the  system  of  equations  (4.4.6)  reduces 


to  the  ubiquitous  laminar  boundary  layer  equations.  The  Prandtl 


hypothesis  holds  and  after  appropriate  nondimensionalization  the 


equations  read 
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(4.4.7a) 
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where  the  streamwise  pressure  gradient  has  been  omitted  because  we 


are  studying  the  flow  on  a  flat  plate  with  zero  pressure  gradient. 


The  primed  variables,  of  course,  correspond  to  the  appropriately 


scaled  boundary  layer  variables,  defined  by  y  =  /Re  y  ; 

Li 


V  =  /Re  V,  where  the  Reynolds  number  is  based  on  the  physical 

^  P„U  L 

properties  of  the  carrier  fluid,  i.e..  Re  =  -  and  L  is  some 


length  scale. 


4.5  The  Turbulent  Regime 


As  is  the  case  in  all  boundary  layer  flow  calculations. 


transition  is  a  crucial  step  of  the  computations.  Its  prediction. 


when  not  known  experimentally,  is  very  difficult, 


Transition  depends  mainly  on  the  free  stream  turbulence  and 


on  the  pressure  gradient.  Since  we  are  doing  computations  on  a 


flat  plate  flow,  which  is  usually  the  first-case  problem  in  fluid 


mechanics,  the  mechanisms  and  position  of  the  onset  of  boundary 
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( 4 . 5 . 3d  ) 


The  attentive  reader  will  have  noticed  that  in  order  to 
obtain  the  equations  above,  disturbances  of  the  void  fraction 
distribution  were  not  considered.  This  avoids  one  of  the  major 
difficulties  in  this  study,  i.e.,  the  simulation  of  the  inter¬ 
action  of  buoyancy  with  the  turbulent  shear  stress. 

If  one  now  performs  an  order  of  magnitude  analysis  of  the 
above  equations  for  a  thin  shear  layer,  i.e.,  in  the  limit  of  high 
Reynolds  number  one  finds  that  the  y-momentum  equation  reduces 
to 
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and  recognizing  the  vanishing  boundary 
y-velocity  fluctuation  and  the  void 


fraction,  we  find 
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If  we  differentiate  this  equation  with  respect  to  x  and  use 
the  fluctuating  counterpart  of  Eq .  (4.5.3b),  we  find  that 
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(4.5.6) 


since  P  =  constant  for  a  flat  plate. 

oo 

It  is  important  to  note  that  the  last  term  of  the  equation 


above  captures,  at  least  in  part,  the  localized  nature  of  the  drag 
reduction.  It  is  clear  that  the  buoyancy  force  field  will  introduce 
a  source  term  in  the  streamwise  momentum  equation.  Furthermore, 


the  importance  of  this  source  term  is  confined  to  the  region  of  the 
flow  where  the  bubbles  are  deflected  by  the  freestream  flow. 

Substitution  of  (4.5.6)  into  (4.5.3c)  eliminates  the  pressure 
term  and  gives  for  the  x  momentum  equation. 
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If  we  estimace  the  order  of  mai^nitude  of  the  various  terms  in 
Eq .  (4.5.7),  we  find  that  in  the  thin  shear  layer,  high  Reynolds 
number  limit  the  x-momentum  equation  reduces  to 
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In  order  to  model  the  turbulent  shear  stress,  we  use  the 
usual  Boussinesq  type  approximation,  i.e.. 


(4.5.9) 


where  e  is  the  eddy  viscosity. 

In  the  inner  region  of  the  boundary  layer,  the  eddy  viscosity 
is  modeled  by  Prandtl's  mixing  length  hypothesis 


(4.5.10) 


where  I,  Prandtl's  mixing  length,  is  computed  from  a  modified 
Van  Driest  expression 

+ 

i  =  khy"^  {1  -  e~^  ,  (4.5.11) 

with  von  Karman's  constant  k  =  0.41. 
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In  this  expression,  followinj;  Launder,  et  al.  [28],  we  use 


A  =  fi/r  (l  +  5.9  v  j  , 


(4.5.12) 


an  expression  that  has  displayed  a  better  behavior  for  flows  with 
boundary  layer  injection.  Through  a  convenient  formulation  of  the 
function  h  we  attempt  to  capture  the  local  nature  of  the  skin 
friction  reduction. 

As  pointed  out  in  Section  4.3,  the  phenomenon  of  microbubble 
boundary  layer  drag  reduction  can  be  consistently  explained  if  one 
thinks  of  it  in  terms  of  momentum  interchanges  between  the  bubbles 
and  the  turbulent  structure  of  the  mean  flow  as  well  as  bubble  flow 
entrainment . 

Under  the  limitations  of  the  description  of  the  flow  adopted 
here,  we  are  led  to  the  conclusion  that  the  mixing-length 
correction  function  h  should  depend  on  the  streamwise  gradient  of 
3v/3y  as  well  as  on  some  nondimensional  parameter  representative 
of  the  average  bubble  radius,  such  as  T^,  i.e.. 


1  +  u  ,2 

h  - - ^  g  T  ' 

1  -  X  ^'3x3y’  a- 


(4.5.13) 


where  the  dependency  on  the  void  fraction  results  from  the  non- 
dimensionalization  of  the  eddy  viscosity  with  the  local  nroperties 
of  the  "fluid". 

The  form  of  the  functional  dependence  of  h  on  the  cross-stream 
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velocity  gradient,  suggested  above,  is  supported  by  the  mathematical 


formulation  of  the  buoyancy  effects. 


It  is  clear  if  the  present  direction  is  adopted  then  the  form 


of  h  requires  further  intensive  study,  the  first  hurdle  of 
which  must  be  the  acquisition  of  a  more  comorehensive  set  of 
experimental  data  on  the  properties  of  the  mean  flow  at  injection. 
For  this  reason,  we  are  forced  to  disregard  g,  even  though  the 
previous  reasoning  suggests  that  g  affords  the  most  practical  way 
of  capturing  the  localized  nature  of  the  drag  reduction. 

In  the  outer  or  wake  region  of  the  boundary  layer,  the  eddy 
viscosity  model  that  we  adopt  is  due  to  Clauser  [29], 

=  aU  6*  ,  (4.5.14) 

0  00 

where  a  =  0.0168  and  6*  is  the  boundary  layer  displacement 
thickness . 

With  expression  (4.5.9)  the  momentum  equation  (4.5.8)  becomes 
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(4.5.15) 


which  along  with  Eq .  (4.5.3a)  and  Eq.  (4.5.3b)  constitutes  the 
system  to  be  integrated. 

The  numerical  scheme  used  to  integrate  the  equations  above  is 
Keller's  Box  Method  [27] .  The  transformation  of  the  above  equations 
by  the  Mangler-Levy-Lees  equations  is  given  in  Appendix  D.  The 
boundary  layer  code  based  on  the  transformed  equations  is  given  in 
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4 . 6  The  Boundary  Condicions 

We  know  from  N.  K.  Madavan's  exDerlmencs  (8]  that  in  the  reeior 
of  bubble  injection  the  flow  is  completely  turbulent.  From  the 
boundary  layer  measurements  in  the  test  section  they  estimate  the 
position  of  a  virtual  flat  plate  leading  edge  that  would  have 
generated  the  boundary  layer  profile  observed.  It  is  based  on  the 
location  of  this  virtual  flat  plate  leading  edge  that  we  decide  the 
location  of  bubble  injection  in  the  code. 

Since  we  are  simulating  the  bubbles  by  a  variable  density  and 
viscosity  "fluid",  we  must  establish  the  boundary  conditions  not 
only  for  the  velocity  field  but  also  for  the  void  fraction.  So, 
if  X  >  0  is  the  streamwise  independent  variable,  we  set 


X(x,0)  =  CX 


ta,b]  ’ 


(4.6.1a) 


where  [a,b]  is  the  interval  of  Injection.  The  estimation  of  the 
constant  C  is  based  on  N.  K.  Madavan's  experiments.  The  laminar 
regime  boundary  condition  for  the  void  fraction  is  given  bv 


X(0,y)  =  0  , 


(4.6.1b) 


since  the  bubbles  are  injected  in  the  turbulent  flow  regime. 

The  velocity  field  boundary  conditions  at  the  plate  are 
determined  by  a  statement  of  conservation  of  mass.  At  the 
impervious  part  of  the  wall  we  have 


V  =  U  =  0  . 


(4.6.2a) 
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In  the  interval  of  injection,  we  impose 

U  =  0  ,  (4.6.2b) 

w 


p  Q  =  V  (b  -  a) [p  f 1  -  X)  +  p^xl 
S  g  w  I-  £  g  J 

or 

g  g 

p^(b  -  a  )  ( 1  -  X) 


(4.6.2c) 


where  Qg  is  Che  volumetric  flowrate  of  the  gas  being  injected. 


4.7  Results 

In  order  to  convince  the  reader  that  the  code  presented  in 
Appendix  E  is  reliable,  we  computed  several  characteristic 
parameters  of  the  boundary  layer  in  the  absence  of  injection. 

An  inspection  of  Figures  4.2,  4.3,  4.4  and  4.5  shows  that 
in  the  absence  of  injection  the  code  generates  results  that 
are  in  agreement  with  widely  known  experimental  data  and 
empirical  relationships  for  the  flat  plate  boundary  layer 
[26,30]. 

Despite  this,  Che  code  is  not  in  a  stage  where  it  can  predict 
void  fraction  profiles  with  high  wall  resolution.  For  a  20  to  30 
grid  point  boundary  layer  resolution,  the  code  predicts  void 
fraction  profiles  that  have  a  maximum  at  y/6  =0.1.  This  fact  is 
very  encouraging  since  this  is  exactly  the  position  on  the 
buoyancy  layer  where  experiments  seem  to  suggest  the  highest 


bubble  concentration. 
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These  results  are  however  not  satisfactory  because  of  the  lack 


of  resolution  of  the  wall  gradients.  The  fact  chat  the  bubble 


injection  occurs  in  Che  turbulent  regime  demands  the  resolution  of 


Che  wall  gradients.  This  can  be  done  rather  easily  in  the  code 


since  the  positions  of  the  first  grid  point  and  the  boundary  layer 


edge  can  be  chosen  freely,  but  judiciously. 


With  this  resolution,  the  coupling  of  the  momentum  and  trans¬ 


port  equations  becomes  very  delicate  because  of  the  high  streamwise 


near  wall  void  fraction  gradients.  Several  unsuccessful  attempts 


were  made  to  correct  this  situation.  A  refinement  of  the  streamwise 


grid  at  injection  as  well  as  a  smoother  injection  void  fraction 


profile  were  unable  to  correct  the  problem. 


The  resolution  of  the  wall  gradients  is  essential  because 


the  form  of  the  transformed  momentum  equation  suggests  that  the 


buoyancy  effects  are  important  only  in  the  inner  region  of  the 


boundary  layer  where  the  turbulent  shear  stresses  are  either 


negligible  or  of  the  same  order  of  magnitude  as  the  viscous 


stresses . 


It  was  observed  that  for  a  rather  modest  wall  void  fraction. 


the  code  before  the  computations  collapsed,  predicts  a  skin 


friction  reduction  of  the  order  of  60%.  This  fact  is  not  only 


encouraging  but  also  points  to  the  possible  solution  of  the 


problems  that  the  model  has  presented. 
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From  Che  beginning  we  assumed  that  the  boundary  layer 


hypotheses  were  valid  throughout  the  boundary  layer  with  micro¬ 


bubbles.  The  large  streamwise  velocity  gradients  at  injection 


suggests  that  the  boundary  layer  hypotheses  are  not  valid  close  to 


the  locus  of  injection. 


These  conclusions  indicate  that  the  order  of  magnitude  analysis 


of  the  momentum  equations  should  be  reviewed. 


The  code  as  it  is  now  should  be  able  to  handle  Che  computations 


of  Che  flowfield  away  from  Che  locus  of  injection,  however,  a 


different  version  of  the  momentum  equations,  as  indicated  above. 


should  be  solved  during  injection. 


This  will  certainly  add  to  the  complexity  of  the  computations 


because  the  momentum  equations  at  injection  will  probably  not  be 


parabolic  anymore  but  elliptic. 
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CHAPTER  V 


DISCUSSION  AND  CONCLUSIONS 

5. 1  The  Stability 

As  was  discussed  in  Chapters  II  and  III,  the  simple  arrangement 
of  hollow  vortices  is  unconditionably  unstable.  It  was  shown, 
through  an  analysis  similar  to  von  (Carman's,  that  the  instability 
is  independent  from  either  the  size  of  the  cylinaers,  as  reflected 
by  the  dipole  strength,  or  the  circulation  of  the  vortices. 

The  evolution  of  a  point  hollow  vortex  sheet  was  included  for 
comparison  purposes.  Even  though  the  time  history  that  we  obtained 
for  this  case  is  different  from  that  presented  by  Chow  [11],  the 
onset  of  the  rollup  is  found  to  have  the  same  pattern.  Indeed,  from 
the  evolution  of  the  hollow  vortex  sheet,  one  concludes  that  the 

consideration  of  the  dipole  system  is  a  destabilizing  factor  since 

the  onset  of  the  rollup  comes  earlier  than  for  the  corresponding 
point  vortex  sheet. 

From  the  point  of  view  of  inviscid  flows,  the  apparent 
stability  of  the  hollow  vortex  sheet  is  associated  with  its 

proximity  to  the  wall.  Indeed,  as  is  clear  from  Figure  2.3,  the 

stability  domain  shrinks  exponentially  away  from  the  wall  and, 
as  expected,  the  condition  of  unconditional  instability  is 
recovered  several  disturbance  wavelengths  away  from  the  wall. 

The  results  of  the  small  disturbance  analysis  of  the 
symmetric  hollow  vortex  sheet  presented  in  Chapter  II  could  not  be 
corroborated  by  the  computation  of  its  evolution  because  the  CPU 
time  required  renders  the  attempt  impractical. 
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For  the  sake  of  completeness,  the  evolution  of  a  sinusoidal 
wave  of  point  vortices  in  the  presence  of  a  wall  was  also  computed. 
Supporting  von  Karman's  work,  the  arrangement  is  always  unstable. 
Indeed,  the  pattern  of  evolution,  vide  Figure  3.3,  is  suggestive  of 
the  large  scale  turbulence  structures  that  develop  in  a  turbulent 
boundary  layer  after  the  breakdown  of  the  Tollmien-Schlichting 
waves . 

It's  conceivable,  and  our  analysis  does  not  address  the  issue, 
that  viscous  effects  also  play  a  role  on  the  stability  of  the 
bubbles,  which  we  simulated  by  cylinders  due  to  the  two-dimensional 
nature  of  the  boundary  layer.  A  deeper  study  of  this  aspect  of  the 
the  problem  is  necessary. 

Attempts  to  correlate,  in  a  very  simple  minded  way  [31],  the 
circulation  of  the  hollow  vortices  to  a  velocity  gradient 
representative  of  their  position  on  the  boundary  layer  lead  to 
meaningless  conclusions.  However,  the  fact  that  the  bubbles  are 
conveyed  by  the  mean  flow,  even  though  not  all  at  the  same  speed, 
renders  the  analytic  treatment  of  the  problem  of  shear  flow 
through  an  array  of  spheres  conceivable. 

The  problem  of  low  Reynolds  number  flow  through  periodic 
arrays  of  spheres  and  cylinders  in  several  oacking  configurations, 
has  been  studied  by  other  authors.  H.  Hasimoto  [32]  studied  the 
problem  of  flow  past  special  cases  of  cubic  lattices.  More 
recently  [33,3<4],  his  work  has  been  expanded,  in  an  equally 
elegant  way,  to  encompass  a  wider  range  of  void  fractions.  The 
latter  authors  derive  an  asymptotic  solution  for  the  high  void 
fraction  flow  regime.  They  provide  the  connection  between  their 
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solution  and  Hasimoto's,  for  the  intermediate  void  fractions, 
through  the  numerical  integration  of  the  resulting  equations 
using  spectral  methods. 

The  extension  of  this  analysis  to  our  problem  must  resolve  the 
difficulty  introduced  by  the  motion  of  the  spheres  convected  by  the 
shear  flow,  i.e.,  the  unsteadiness.  The  unsteadiness  translates 
itself  into  a  time  dependency  of  the  solution  domain.  This  is 
where  the  mathematical  difficulties  arise. 


5.2  Drag  Reduction 

If  not  from  the  title  of  this  dissertation  certainly  from 
Chapter  IV,  the  reader  must  have  realized  the  monumental  task  that 
it  will  be  to  describe  this  flow  even  from  an  engineering  point  of 
view. 

Ahead  we  try  to  point  out  further  some  more  of  the 
difficulties  and  suggest,  in  broad  lines,  some  research  that 
can  lead  to  a  better  understanding  of  the  problem. 

The  fact  that  the  bubbles  are  very  small  is  very  suggestive 
of  a  phenomenological  approach  to  this  problem.  On  the  other 
hand,  the  high  packing  of  the  bubbles  may  support  arguments  in 
the  opposite  direction.  Whichever  direction  one  choses,  the 
challenge  is  enormous. 

If  one  decides  to  formulate  appropriate  constitutive  relations 
for  a  bubbly  flow  [22],  then  the  present  experimental  techniques 
must  be  improved  to  allow  the  measurement  of  mean  flow  properties 
in  regions  of  high  void  fraction.  To  the  author's  knowledge,  very 
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little  effort  has  been  conducted  in  this  direction,  at  least  in 
connection  with  this  problem.  Nevertheless,  exciting  progress 
seems  to  have  been  achieved  recently  by  R.  van  der  Welle  [35].  He 
reports  the  development  of  ingenious  measurement  techniques  for 
bubbly  flows  with  high  void  fraction  that  can  conceivably  be 
extended  to  the  study  of  the  present  flow. 

It  is  important  to  notice  that  any  such  constitutive  relation 
will  be  dependent  not  only  on  local  properties  of  the  fluid  but 
also  on  local  properties  of  the  flow  [36,37,38].  The  latter  effect 
is  dictated  by  the  dynamic  interaction  between  the  bubbles  and  the 
carrier  flow. 

Another  very  important  aspect  of  an  experimental  investigation 
must  be  the  study  of  the  interaction  of  the  turbulent  structure 
and  the  carrier  flow  with  the  dynamics  of  bubbles.  This  complex 
phenomenon  seems  to  be  very  important  at  least  in  the  region  of 
injection  of  the  bubbles.  Indeed,  the  Soviet  researchers  [2,4] 
postulate  that  this  is  the  most  important  mechanism  of  the  drag 
reduction.  If  this  is  indeed  the  whole  picture,  it  is  not  clear. 
The  idea  is  supported  however,  by  the  fact  that  the  drag  reduction 
persists  only  for  a  few  boundary  layer  thicknesses  downstream 
from  the  place  of  injection.  The  way  in  which  this  may  be 
happening  was  presented  in  Section  4.3. 

One  is  faced  with  a  major  comnutational  and  modeling  effort, 
specially  in  providing  appropriate  closure  models  for  the  resulting 
turbulent  flow  equations.  Attention  must  be  dedicated  not  only  to 
the  modeling  of  the  buoyancy  effects  [19,39]  but  also  to  the  shear 
stresses  themselves.  It  is  not  clear,  for  example,  chat  the 


normal  stresses  in  the  streamwise  and  cross-stream  directions  are  of 
the  same  order  of  magnitude  as  is  usually  assumed,  and  is  experi¬ 
mentally  confirmed  in  the  flat  plate  turbulent  boundary  layer. 
Another  aspect  is  that,  since  buoyancy  and  the  void  fraction  play 
such  an  important  role  on  this  flow,  it  is  probably  more  appropriate 
in  the  averaging  of  the  Reynolds  number  to  use  a  mass  weighed 
ensemble  average  [39]. 

Also,  a  mixing  length  approach  in  the  modeling  of  the  shear 
stresses  may  prove  inadequate.  The  empirical  relations  used  in 
this  approach  depend  on  particular  characteristics  of  the 
boundary  layer  which  have  not  yet  been  confirmed  true  for  the 
boundary  layer  with  microbubbles.  This  suggests  interesting 
experimental  research  that  could  make  use  of  the  techniques 
mentioned  in  [35].  The  strong  interaction  of  the  bubbles  with 
the  turbulence  structure  may  preclude  any  such  approach  and  force 
the  use  of  second  order  closure  models  [39].  It  is  not  clear, 
for  example,  that  the  interraittency  characteristics  exhibited  by 
boundary  layers  with  and  without  bubbles  is  the  same.  Indeed, 
the  reasoning  presented  before  on  the  mechanisms  of  the  drag 
reduction  suggests  that  this  is  not  the  case.  If  experiments 
ever  substantiate  this,  then  an  appropriate  formulation  for  the 
eddy  viscosity  in  the  outer  region  of  the  boundary  layer  is 
necessary. 

Either  in  the  approach  adopted  here  or  In  a  complete  two 
phase  flow  description,  the  researcher  must  be  concerned  about 
the  appropriate  formulation  of  the  boundary  conditions.  The 
carrier  flow  boundary  conditions  are  straightforward,  however. 


it  is  clear  that  the  wall  boundary  condition  for  the  void  fraction 
is  itself  a  function  of  the  injection  flow  rate.  Consequently, 
even  though  the  approach  adopted  here  is  physically  meaningful, 
it  is  rather  restricted  and  must  be  improved.  Once  again, 
measurement  techniques  must  be  perfected  in  order  to  allow  near 
wall  void  fraction  measurements  and  determine  its  correlation  with 
the  injection  flowrate. 

For  the  sake  of  completeness  and  adding  to  the  complexities 
already  mentioned,  the  author  must  point  out  that  recent  results 
from  ongoing  research  in  this  problem  [40]  suggests  that  the 
bubbles  are  bigger  than  originally  thought  and  the  coalescence 
is  a  fact  specially  downstream  from  injection. 

It  is  clear  from  what  we  presented  as  well  as  past  work 
[2-81  that,  at  least  in  the  near  future,  the  main  effort  towards 
understanding  this  problem  must  be  carried  out  experimentally. 

At  the  current  research  stage  there  are  too  many  gaps  in  relevant 
experimental  data.  A  more  complete  knowledge  of  these  is  needed 
if  the  researcher  is  to  approach  the  modeling  of  this  flow  with  a 
reasonable  degree  of  confidence. 
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APPENDIX  A 


THE  FOURIER  SERIES  EXPANSIONS 


In  order  to  obtain  expressions  (2.3.8),  several  Fourier  series 
expansions  were  obtained: 


coshk.0  =  sinhkTT  +  T  „  sinhkn  cosmr  cosn0  ,  (A.l) 

Trk  ^2,2 

n=i  n  +  k 

"  .2  _  2v 
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A  combination  of  expressions  (A.l)  with  (A. 2)  and  of  (A. 3)  with 
(A. 4)  leads  to  Eq .  (2.3.8d)  and  Eq .  (2.3.8a),  respectively. 

Equation  (2.3.8d)  is  obtained  from  expressions  (A. 4),  (A. 5)  and  (A. 6) 
while  Eq .  (2.3.8c)  stems  from  (A.l),  (A. 7)  and  (A. 8). 


APPENDIX  B 


THE  CONTOUR  INTEGRALS 


To  integrate  Eq .  (3.3.17)  we  use  Cauchy's  Residue  Theorem.  In 
order  to  perform  the  integrations  we  need  the  following  expressions 
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The  integral  /  dz  leads  to  integrals  of  the  type 
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The  integral  /  leads  to  integrals  of  tl  type 
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dW 

The  integral  /  "Y”  dz  leads  to  integrals  of  the  type  (B.3) 
Cm  " 

and  (B.4).  With  these  results  and  taking  into  consideration  the 
appropriate  coefficients,  Eqs.  (3.3.8)  can  be  readily  obtained. 
For  the  image  system,  the  expressions  to  integrate  in 
Eq.  (3.3.17)  are 
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From  these  expressions  it  is  obvious  that  in  Blasius  equations 
only  Eq.  (3.3.18b)  generates  more  terms. 

With  expression  (B.7),  Eq .  (3.3.18b)  produces,  in  addition 
to  (B.5)  and  (B.6),  integrals  of  the  following  type 
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With  these  results  at  hand,  Eq .  (3.4.5)  follows  easilv. 


IMPLICIT  REAL*8(A-H,0-Y) 

IMPLICIT  C0MPLEX*16 ( Z ) 

COMMON/BLKl/J , N , P I 

COMMON/BLK2/ZPOS 

COMMON/BLK3/ZVEL 

COMMON/BLK4/DT 

COMMOK/BLKS/ZPC 

DIMENSION  ZPOS(IOO) ,ZPOS0(50) ,ZPOSl( 50) 
DIMENSION  ZVELdOO) 

PRINT*, 'ENTER  N,M' 

READ ( 5 , * )  N , M 
PRINT* ENTER  AD, DEL' 

READ{5,*)  AD, DEL 

Nl-N+1 

N2«N/2 

N21-N2+1 

DT-l/DFLOAT( 500 ) 

PI=DACOS(-l .DO) 

SI=DCM?LX( 0 .DO , 1 .DO  ) 

DO  10  I»1,N 

X-( I-l )/DFLOAT(N) 

Y=-AD*DSI.N(2*?I*(X)  ) 

ZPOS( I )=DCMPLX(X,Y) 

WRITE(21,*)  X,Y 
10  CONTINUE 

ZVEL ( 1 ) -DCMPLX ( 0 . ODO , 0 . ODO ) 

ZVEL  ( .N21 )  =DCMPLX  (  0  .  ODO  ,  0  .  ODO  ) 

DO  2  0  J»2,.N2 
ZPC»ZPOS( J) 

ZVEL( J)»2I*(DCONJG(ZSCOT(J,ZPOS( J) ) ) )/N 
ZVEL(N+2-J)*«-ZVEL(  J) 

20  CONTINUE 
DO  60  K=1,M 
T-T+DT 

CALL  ORDEROv ZFOSO ) 

CALL  ORDERl ( ZPOSl ) 

DO  30  1=2, N2 

ZFCS( I )=ZFOS0 ( I )  +  (DEL**2 ) *ZP0S1 ( I  ) 
ZP0S(N+2-I )=1 .D0-ZPOS( I ) 

30  CONTINUE 

DO  40  J=2,N2 
ZPC=ZPCS ( J ) 

ZVEL( J)=ZI*{DCONJG{ZSCOT{ J,ZPOS( J) ) ) ) /N 
ZVEL(N+2-J )--ZVEL( J ) 

40  CONTINUE 

IF(MOD( DFLOAT( K) , 4 0 . DO ) . EQ . 0 . DO )  THEN 

L1-21+K/40 

DO  50  1=1, N1 

X-DREAL( ZPOS( I ) ) 

Y-DIMAG{ ZPOS( I ) ) 

WRITE{L1,*)  X,Y 
50  CONTINUE 

WRITE(L1,*)  T 
END  IF 

60  CONTINUE 


ZD4-DT*ZRHS1( ZPC) 

ZPOUT-ZPIN+DCONJG( ZD1+2*ZD2+2*ZD3+ZD4 )/6 . ODO 

RETURN 

END 


SUBROUTINE  RUNTA2 ( DT , ZPIN , ZVIN , ZPOUT ) 

IMPLICIT  REAL*8 ( A-H ,0-Y) 

IMPLICIT  COMPLEX*16 ( Z ) 

COMMON/BLKS/ZPC 
ZDV1=DT*ZRHS2 ( ZPIN, ZVIN) 

ZDP1=DT*ZRHS3 ( ZVIN) 

ZPC-ZPIN+DCONJG( ZDPl )/2 . ODO 
ZVC=ZVIN+DCONJG( ZDVl )/2 . ODO 
ZDV2-DT*ZRHS2 ( ZPC,ZVC) 

ZDP2=DT*2RHS3 ( ZVC ) 

ZPC=ZPIN+DC0NJG( ZDP2 )/2 . ODO 
ZVC=ZVIN+DCONJG( ZDV2 )/2 . ODO 
ZDV3=DT*ZRHS2 ( ZPC , ZVC ) 

ZDP3=DT*ZRHS3 ( ZVC ) 

ZPC=ZPIN+DC0NJG{ ZDP3 ) 

ZVC-ZVIN+DCONJG( ZDV3 ) 

ZDV4*DT*ZRHS2( ZPC,ZVC) 

ZDP4-DT*ZRHS3 ( ZVC ) 

ZPOUT=ZPIN+DCONJG( ZDP1+2*ZDP2+2*ZDP3+ZDP4 ) /6 . ODO 

RETURN 

END 


FUNCTION  ZRHSl(ZPC) 
IMPLICIT  REAL*8( A-H,C-Y) 
IMPLICIT  COMPLEX*16 ( Z ) 
COMMON/BLKl/J ,N, PI 
ZI-DCMPLX( 0 .DO , 1 .DO ) 
ZRHSl='-ZI*ZSCOT(  J,ZPC)/N 
RETURN 
END 


FUNCTION  ZRHS2 ( ZPC, ZVC) 
IMPLICIT  REAL*8( A-H,0-Y) 
IMPLICIT  C0MPLEX*16 ( Z ) 
COMMON/BLKl/J , N , PI 
COMMON/BLK2/ZPOS 
CCMMC^?  /BLK  3  ,''ZVEL 
DIMENSION  ZPOS( 100 ) ,ZVEL( 100) 
ZI-DCMPLX( 0.D0,1.D0) 
ZS1=ZSC0T( J,ZPC) 
ZS2-DCMPLX(0.D0,0.D0) 


STOP 

END 


SUBROUTINE  ORDERO ( ZPOSO ) 

IMPLICIT  REAL*8(A-H,0-Y) 

IMPLICIT  COMPLEX*16{Z) 

COMMON/BLKl/J,N, PI 

COMMON/BLK2/ZPOS 

COMMON/BLK4/DT 

COMMON/BLK5/ZPC 

DIMENSION  ZPOS0(50) ,ZPOS{100) 

N2=N/2 

DO  10  J=2,N2 
ZPC=ZPOS( J) 

CALL  RUNTAl(DT,ZPOS( J) ,ZPOS0( J) ) 
CONTINUE 
RETURN 
END 


SUBROUTINE  ORDERl ( ZPOSl ) 

IMPLICIT  REAL*8 (A-H ,0-Y) 

IMPLICIT  COMPLEX*16 ( Z ) 

COMMON/BLKl/J,N,PI 

COMMON/BLK2/ZPOS 

COMMON/BLK3/ZVEL 

COMMON/BLK4/DT 

COMMON/BLK5/ZPC 

DIMENSION  ZPOS(IOO) , ZPOSl (50) ,ZVEL(100) 
N2-N/2 

DO  10  J=2,N2 
ZPC=ZPOS ( J ) 

CALL  RUNTA2(DT,ZPOS( J)  ,ZVEL( J)  ,ZPOSl( J)  ) 
CONTINUE 
RETURN 
END 


SUBROUTINE  RUNTAl ( DT , ZPIN , Z POUT ) 
IMPLICIT  REAL*8{A-H,0-Y) 

IMPLICIT  COMPLEX*16(Z) 
COMMON/BLK5/ZPC 
ZD1-DT*ZRHS1 ( ZPIN) 

2?c-3r;;;hDccr:jG{  zdi  .  odo 
ZD2-DT*ZRHS1 ( ZPC ) 
ZPC-ZPIN+DCONJG( ZD2 )/2 . ODO 
ZD3=-DT*ZRHS1(ZPC) 
ZPC=ZPIN+DCONJG( ZD3 ) 


ZS3=-DCMPLX(  0  .DO  ,  0  .DO  ) 

DO  10  1=1, N 

IF( I .EQ. J)  GO  TO  10 
Z-PI* ( ZPC-ZPOS ( I  )  ) 

ZS2  =  ZS2+(ZS1-ZSC0T( I,ZPOS( I  )  )  )/(  (ZDSIN(Z)  )**2) 
ZS3=ZS3+( ZVC-ZVEL( I ) )/( ( ZDSIN( Z ) ) **2 ) 

10  CONTINUE 

ZS2=4*ZI*ZS2/( N**2 ) 

ZS3=-2* ( ZI*ZS3+2*DIKAG( ZS3 ) )/N 
ZRHS2  =  PI* ( ZS2  +  ZS3  ) 

RETURN 

END 


FUNCTION  ZRHS3(ZVC) 
IMPLICIT  REAL*8 (A-H,0-Y) 
IMPLICIT  C0MPLEX*16 (Z  ) 
ZRHS3=DCONJG( ZVC ) 

RETURN 

END 


FUNCTION  ZSCOT(L,ZP) 
IMPLICIT  REAL*8(A-H,0-Y) 
IMPLICIT  COMPLEX*16(2) 
COMMON/BLKl/J,N, PI 
COMMON/BLK2/ZPOS 
COMMON/BLK5/ZPC 
DIMENSION  ZPOS(IOO) 
ZS-DCMPLX( 0 .DO , 0 .DO  ) 

DO  10  1  =  1, N 

IF(I.EQ.L)  GO  TO  10 
IF(I.EQ.J)  THEN 
ZPl=Z?C 
ELSE 

ZPl=ZPOS { I ) 

ENDIF 

Z  =  PI* ( ZP-ZPl  ) 

ZS-ZS  +  ZDCOS ( Z )/ZDSIN( Z  ) 
10  CONTINUE 
ZSCOT=ZS 
RETURN 
END 


FUNCTION  ZDSIN(Z) 
IMPLICIT  REAL*8 ( A-H,0-Y) 
IMPLICIT  C0MPLEX*16 ( Z  ) 
ZI-DCMPLX{ 0 . DO , 1 .DO  ) 


XR=-DREAL(  Z  ) 

XI-DIMAG( Z ) 

X-DSIN(XR) *DCOSH(XI  ) 

Y-DCOS(XR) *DSINH(XI  ) 

ZDSIN-X+ZI*y 

RETURN 

END 


FUNCTION  ZDCOS(Z) 
IMPLICIT  REAL*8 (A-H,0- 
IMPLICIT  CCMPLEX*16 ( Z  ) 
ZI-DCMPLX( 0 . DO , 1 .DO ) 
XR-DREAL( Z ) 

XI-DIMAG( Z ) 

X-DCOS{XR)*DCOSH{XI ) 

Y-DSIN(XR) *DSINH(XI ) 

ZDCOS-X-ZI*Y 

RETURN 

END 


97 


APPENDIX  D 


THE  MANGLER-LEVEL-LEES  TRANSFORMATION 


In  Chapter  IV,  the  system  of  equations  to  be  integrated  is 


3u  3v  „ 
—  +  —  =  0 
3x  3y 


(D. la) 


(Xu)  +  (Xv)  =  0 
3x  3y 


(D. lb) 


X,  r  3u  1  /.  3 _  r  *  3u  1 

(1  -  x)[u  —  +  V  — ]  .  (1  -  X)  —  [VE 


V  ^ 


(D.lc) 


where  e*  =  e/  v. 

If  we  define  a  stream  function  ip  according  to 


3ip  3ip 

^  “  3y  ’  ~  3x  ’ 


(D.2) 


then  the  continuity  equation  is  identically  satisfied. 

With  a  Mangier -Levy-Lees  type  of  transformation,  defined  by 


C(x)  =  U_^x 

ref 


(D.3a) 


n(x,y) 


(D.3b) 


ip(C,n)  =  /2c  f  (C.n)  , 


(D.3c) 


I 

'Wi 


the  two  remaining  equations  become 


2,  f.||-(f  =0  , 


( D. 4a ) 


-^1-1(1  *  ^J«"l  +  <i  -  «  - 

li  8n  ^  sp-'  J  u 

ref  ref 


^  "sp^  3  re*f" 


d  rEL  7 

ll-ril  +  ”  1) 


^25 


t‘  *  25  Ij-  nijK/  XdO 


(1  -  X)25[f ^f '  -  f " 


(D.4b) 


The  corresponding  boundary  conditions  at  the  wall  become 


f  (C,0)  =  0  , 


(D.5a) 


/25  f(?,0)  =  -— /  V  (Oc“C  , 
ref 


(D.5b) 


while  the  farfield  boundary  condition  is 


f(?,n)-^lasri+”  . 


(D.5c) 


The  volume  fraction  boundary  conditions  remain  as  expressed 
in  Eq.  (4.6.1). 

In  order  to  apply  Keller's  Box  Method,  we  rewrite  the  system 
of  equations  as 


F  =  f 


G  =  f'  =  f  '  ' 


(D.6a) 


(D.6b) 


2CFX^  -  (f  +  2Cf^)Y  =  0  , 


(D.6c) 


— ^  [(1  +  U  ]g]  +  (1  -  X)  {-^ - 

u  sp^  •>  u. 

ref  "^ref 


J 


-L/2J  [1  ^  25!^.  „i.J(J  Xdc) 

Pj^U*  n 


=  (1  -  X)2c[FF^  -  f^G]  , 


(D.6d) 


vd.th  the  boundary  conditions 


F  =  0 


at  n  =  0  , 


F  -••  1 


X  =  CX„ 


-  -L- 

/U  ^2  “  ^1 


as  n  ■*•  ®  , 


at  n  =  0  , 


at  n  =  0  , 


(D.6e) 


where  D  is  the  interval,  in  the  transformed  space  (C»n), 
corresponding  to  [a,b]. 

For  the  sake  of  completeness  we  give  now,  in  transformed 
variables,  the  expressions  used  for  the  eddy  viscosity. 

In  the  inner  region  we  have 


—  /2C  f ' 'n^[l 


+  0  1/2 

Z_  (Jf. 

^  P  T  ^ 

A  w 


(D.7a) 


where 


(D.7b) 
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1  -  X  U  , .  . 

w  f  . 
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In  the  outer  region  we  used 


e„  =  a  ^  f  -  ( I  -  X)f ' ]dn 
°  0  ^ 


(D.7c) 


(D.7d) 


Finite  differencing  the  system  (D.6)  according  to  Keller's 
Method  [27],  leads  to  a  block  tridiagonal  system,  with  3x3  block 
matrices,  for  the  parabolic  momentum  equation.  The  volume  fraction 
equation,  a  hyperbolic  equation,  is  integrated  separately,  in  each 
cycle,  by  a  different  solver.  The  differencing  scheme  chosen  for 
this  equation  is  second  order  accurate. 

Since  the  numerical  methods  are  not  the  object  of  this 
dissertation,  we  decline  from  including  the  finite  difference 
equations  that  lead  to  the  code  presented  in  Appendix  E. 


APPENDIX  E 


THE  BOUNDARY  LAYER  CODE 

The  code  presented  in  this  appendix  was  designed  for  the  flat 
plate  boundary  layer.  It  corresponds  to  Che  numerical  integration 
of  the  system  (D.6). 

As  it  is,  the  code  cannot  be  used  for  the  computation  of  the 
flow  around  an  axis3nnmecric  body.  However,  the  changes  to  be 
introduced  are  rather  straightforward.  In  this  case,  the  Reynolds 
equations  expressed  in  body  fitted  coordinates  should  be  trans¬ 
formed  through  the  Stokes  stream  function  and  Che  axisymmetric 
Manger-Levy-Lees  transformation  Co  obtain  the  self-starting  system 
of  equations  [27].  The  system  obtained  in  the  previous  appendix 
must  be  the  special  case  of  zero  pressure  gradient.  Some  entries 
of  the  matrices  of  the  block  cridiagonal  system  will  be  different, 
however,  the  numerical  solver  is  exactly  the  same  and  needs  no 
modifications.  It  should  be  noted  that  the  body  shape  will  also 
affect  both  equations  by  the  definition  of  the  Stokes  stream 
function  and  the  Mangler-Levy-Lees  transformation. 

Even  though  confident  about  the  code  presented,  the  author 
wants  to  issue  a  word  of  caution.  The  ability  of  the  code  to 
resolve  high  void  fraction  gradients  is  limited. 

The  code  uses  second  order  differencing  everywhere  and  is 
designed  for  a  variable  grid  in  either  the  screamwise  or  the 
cross-stream  directions. 


Because  the  void  fraction  transport  eauation  introduces  in 
the  flow  strong  streamwise  gradients,  there  exists  the  need  of  a 
grid  generator  in  the  streamwise  direction.  This  addition  can 
be  appended  to  the  code  without  any  modifications. 


non 


*  MAIN  PROGRAM  * 

★  * 


IMPLICIT  REAL*8 ( A-H , 0-Z ) 

DIMENSION  A(3,3,101),B(3,3,101),C(3,3,101),R(3,101) 

DIMENSION  AL( 3,3, 101 )  ,BET( 3,3,101 ) ,DETA( 101 ) , BUOY ( 101  ) 

DIMENSION  Y(3,101),Z(3,101), ITER{1000) ,ALFA( 1000 ) ,V( 3,1 ) 

DIMENSION  SFI(101),BFI(101),GI(101),ZI(1000), XI 1(101) 

DIMENSION  SFI1(101),BFI1(101),GI 1(101) ,WKAREA( 10 ) 

DIMENSION  AHAT( 101 )  ,BHAT( 101 ) ,CHAT( 101 )  ,DHAT( 1101  ) 

DIMENSION  P(101),Q(101),S(101) ,REX( 1000) ,YDEL( 101) 

DIMENSION  SFI2(101), BFI 2(101), GI 2(1 0-1), CF (1000), XI (101) 

DIMENSION  THMO( 1000 ) ,XDEL( 500 ) ,QBL( 500 ) , XI 2 ( 1 0 1 ) , DU ( 3 , 3 ) 

DIMENSION  UPLUS( 101 )  , YPLUS( 101 ) , SUM ( 101 )  , SUMl( 101  ) 

DIMENSION  SUM2( 101 ) ,BUOYl( 51) , SHAPE ( 1000 ) , ETA( 101 ) 

DIMENSION  PXL(lOl) ,PXL1(101) ,THBL(1000) ,SFI3(101) 

COMMON/BLKl/DETA, CE,N 

COMMON/BLK2/ALFA, IS 

COMMON/3LK3/SFI ,BFI  ,GI 

COMMON/BLK4/SFI1 ,BFI1 ,GI1 

COMMON/BLK5/PRTR , PRTRl , PRST, PRSTl 

COMMON/BLK6/XI , XI 1 

COMMON/BLK7/Z I , ETA 

COMMON/BLK8/PXL , PXLl 

COMMON/BLK9/ISINJI , EREF , VI SCL , VI SCT 

COMMON/BLKl 0/GRAV , ROUL , UE 

COMMON/BLKll/SUM, SUMl ,SUM2 

COMMON/BLK12/YPLUS ,USTAR 

COMMON/BLKl 3/SFI 2 , XI 2 

COMMON/BLKl 4/BUOY , BUOYl 

STRFAC-DLOG( 10 . ODO ) 

EPS-1 . OD-5 
PI-DACOS(-l , ODO ) 

VISCL-lOlOD-6 
VISCG-17 , 9D-6 
ROUL-997 . 3 
ROUG-1 . 21 
GRAV-9 .81 
PRINT* , ' ENTER  M' 

READ(5,*)  M 

Specification  of  flow  parameters. 

PRINT* ,' ENTER  INITIAL  POSITION' 

READ(5,*)  PIIC 

PRINT* ,' ENTER  FINAL  POSITION' 

READ(5,*)  PFIC 
PRINT* , ' ENTER  VFI , QG' 

READ(5,*)  VFI,QG 


PRIimT*  ,  '  ENTER  UE' 

READ(5,*)  UE 

ISINJI-0 

ISINJF-0 

POSEND=0 . ODO 

DZIC-1 . ODO/DFLOAT(M) 

DZI-DZIC 
VISCT=VISCL 
PRTR-0 . ODO 
PRST=1 . ODO 
POS=0 . ODO 
CE=1 . 05D0 
IS--1 

YDEL( 1 )-0 . ODO 
ALFA(  1  )>-0  .  ODO 

Establishment  of  the  cross  stream  grid  which  obeys 
a  geometric  progression  law. 

ETA(1)-0.0D0 
DETA( 1)-O.ODO 
DETA( 2 )-0 .005 
ETA( 2 )*DETA( 2 ) 

K-2 

Determination  of  the  number  of  grid  points  across 
the  boundary  layer. 

1  DETA(K+1)-DETA(K)*CE 

ETA( K+1 )-ETA( K)+DETA( K+1 ) 

ETAE-DETA( 2)*(CE**(K-1)-1.0D0)/(CE-1 .ODO  ) 
IF(ETAE.LE.8.0D0)  then 
N-K 
K»K+1 
GO  TO  1 
ENDIF 
NP1=N+1 

Specification  of  the  starting  parabolic  profile. 


ZI ( 1  )=0  . 
SFI ( 1 )-0 
BFI (1 )=0 
XI ( 1 )-0  . 
XII ( 1 )-0 
GI( l)-2. 
YDEL( 1 )- 
BUOY( 1 )- 
BUOYl ( 1 ) 
SUM( 1 )-0 
SUMl ( 1 )- 
SUM2( 1  )- 
PXL( 1 )-0 
PXLl ( 1 )- 
SFI3( 1 )- 


ODO 
.  ODO 
.  ODO 
ODO 
.  ODO 

ODO/ETAE 
0  .  ODO 
0  .  ODO 
-0 . ODO 
.  ODO 
O.ODO 
0  .  ODO 
.ODO 
0  .  ODO 
0  .  ODO 


Ct 

S 


I 


DO  5  K-2,NPl 

IF(ETA(K) .LE. { 5.0D0) )  IDEL=K 
SFK  K)-(ETA( K) **2 ) * ( 1 . 0D0-ETA( K)/ 

> ( 3*ETAE ) )/ETAE 

BFI { K)-ETA( K) * ( 2 . 0D0-ETA(  K )/ETAE )/ETAE 

GI  ( K)-2*( 1 . 0D0-ETA( K)/ETAE)/ETAE 

XI ( K)-0 . ODO 

XI1(K)=O.ODO 

BUOY( K)=0 . ODO 

BUOY1(K)=O.ODO 

SUM( K)»0 . ODO 

SUM1(K)<»0.0D0 

SUM2(K)-0.0D0 

PXL(K)=0. ODO 

PXL1(K)-0.0D0 

SFI3(K)»0.0D0 

YDEL(K)-BFI(K) 

5  CONTINUE 

10  IF(IS . EQ. ISINJI )  THEN 
DO  15  X-1,NP1 
SFI2(K)=SFI3(K) 

SFI1{K)-SFI2(K) 

SFI(K)«SFI2(K) 

BFI1(K)“BFI2(K) 

BFI(K)=BFI2{K) 

GI1(K)-GI2(K) 

GI ( K)-GI2( K) 

15  CONTINUE 
END  IF 

Computation  of  the  streamwise  Reynolds  number. 
POSl-POS 

POS=POS+DZI/VISCT 
20  IFdS.EQ.l)  THEN 
REX(  1  )=i0,  ODO 
ELSE 

REX( IS )-POS/VISCL 
ENDIF 

Computation  of  the  transition  parameter  that  switches 
on  the  turbulence  model, 

IF(REX(IS) ,LT.4.D5)  THEN 
PRTR-O . ODO 

Computation  of  the  laminar  boundary  layer  edge 
parameters . 

PRCF-DSQRT( 2 . ODO ) *GI { 1 ) 

PRSFI-( SFI ( IDEL+1 )-SFI  (  I DEL)  )/DETA( I  DEL ) 

PRBFI-( BFI ( IDEL  +  1 )-BFI  {  I DEL)  )/DETA( I  DEL) 

SFIID-SFK IDEL)+PRSFI*( 5. ODO- ETA ( I DEL) ) 

BFIID-BFK IDEL)+PRBFI*  (  5 . ODO- ETA ( I DEL )  ) 

PREV— (SFIID-5. 0D0*BFIID)/DSQRT(  2  .  ODO  ) 


ELSE 

IF(REX(IS )  .LE.6 .05)  THEN 
X=REX( IS )/l .05 
PRTRl-PRTR 
X-( X-4 . OOO )/2 . OOO 

PRTR-( 1 . 0D0+OSIN( PI* (X'O . 500 ) ) )/2 . OOO 

PRTR=PRTR**2 

ELSE 

PRTR1=PRTR 
PRTR=1 . OOO 
ENDIF 
ENDIF 

Computation  of  the  boundary  layer  displacement  thic)<ness 
and  of  the  boundary  layer  momentum  thic)cness. 

IF(IS.GT.2)  THEN 

CALL  BLPR(OELSTR,TETSTR) 

THBL( IS-1 )=OSQRT( 2*21 { IS-1 ) ) *OELSTR/( UE*ROUL ) 

THMO( IS-1 )=OSQRT( 2*ZI ( IS-1 ) ) *TETSTR/( UE*ROUL ) 

Computation  of  the  boundary  layer  shape  factor. 

SHAPE  (  IS-1  )=«THBL(  IS-1  )/THMO(  IS-1  ) 

ENDIF 

PRST1=PRST 

P0SC1=P0S1/(UE*R0UL) 

POSC=POS/(UE*ROUL) 

Determination  of  the  final  station  of  injection  talking 
into  account  the  influence  of  the  void  fraction. 

IF((IS.GT.ISINJF) .AND. ( ISINJF.NE. 0 ) )  THEN 
IF(POSC.LE.PFIC)  THEN 
ISINJF=IS 
ETAREF=»ETA(  NPl  ) 

ELSE 

Computation  of  the  position  downstream  from  injection 
in  terms  of  the  boundary  layer  thickness. 

XDEL(  IS-1  )=POSEND/(  ETAREF*DSQRT  (  2*PFI.H  )/ 

>(UE*ROUL) > 

ENDIF 

ENDIF 


Determination  of  the  initial  station  of  injection. 

IF( ( PIIC .GE. POSCl ) .AND. ( PI IC . LE . POSC ) )  THEN 
IF( ISINJI .NE. 0 )  GO  TO  30 
ISINJI-IS 

PIIM-(ZI(IS)+ZI(IS-1))/2.0D0 

POSIN-POSl 

DZIIN-DZI 

VISCTI=VISCT 


'-■vV- 


PRTRI-PRTR 

PRSTI=PRST1 

NIN-N 

C 

C  Storage  of  the  converged  initial  injection  profile 

C  to  restart  the  computations  once  the  final  position 

C  of  injection  is  determined. 

C 

DO  25  K=1,NP1 

SFI 3 ( K ) =SFI2 ( K ) 

SFI2(K)=SFI(K) 

BFI2(K)=BFI(K) 

GI2(K)=GI(K) 

25  CONTINUE 
30  ENDIF 
C 

C  Determination  of  the  final  station  of  injection. 

C 

IF( (PFIC.GE.POSCl) .AND. ( PFIC . LE . POSC ) )  THEN 
IF( ISINJF.NE.O)  GO  TO  35 
ISINJF=IS 
ETAREF-ETA(NP1 ) 

PFIM-ZI ( IS )+DZI/2 . ODO 
IS-ISINJI 
POS=POSIN 
DZIC-DZIIN 
VISCTF=VISCT 
VISCT=VISCTI 
PRTR-PRTRI 
PRST-PRSTI 
N-NIN 
NPl-N+1 
C 

C  Restart  of  the  computations  from  the  initial  position  of 

C  injection.  With  the  final  position  of  injection  in  the 

C  transformed  space  determined,  injection  and  the  transport 

C  equation  for  the  void  fraction  will  be  turned  on. 

C 

GO  TO  10 
35  ENDIF 
C 

C  These  commands  change  the  boundary  conditions  during 

C  injection  according  to  a  statement  of  mass  conservation. 

C 

IF({ZI(IS).GT,PIIM) .AND.(ZI(IS) .LT.PFIM) )  THEN 
PRPOS=-(  PFIM-PIIM)/10  .ODO 

PRINJ-ROUG*QG/( ( PFIM-PI IM ) *DSQRT( 2.0D0*ZI(IS))) 

IF( ZI ( IS ) . LE( PIIM+PRPOS ) )  THEN 
XI ( 1 )-VFI* ( ZI( IS )-PIIM)/PRPOS 
ELSE 

IF(ZI{ IS)  .GE. (PFIM-PRPOS)  )  THEN 

XI ( 1  )-VFI-VFI* ( ZI ( IS )-PFIM  +  PRPOS )/PRPOS 
ELSE 

XI  { 1  )-VFI 
ENDIF 


■/s’' 


nnnn  ooo  non  noon  onno  noonn 


ENDIF 

VOLBC=VOLBC+DZI* ( 1 . 0D0/{ l.ODO-XI(l) )+1.0D0/(1.0DO 
>-XIl ( 1 ) ) )/2 . ODO 

SFIBC=-PRINJ*VOLBC 

ENDIF 

IF((ZI(IS).GE.PFIM) .AND. (PFIM.NE.O.ODO))  THEN 
XI(1)=0.0D0 

SFIBC  =  -VOLBC*ROUG*QG/(  (PFIM-PIIM) *DSQRT{ 2*ZI ( IS  )  )  ) 
ENDIF 

Initialization  of  the  void  fraction  profile  az  the 
initial  station  of  injection  in  order  to  satisfy  the 
zero  wall  void  fraction  gradient  condition. 

IF( IS . EQ. ISINJI )  THEN 

XI(2)=XI{1)*CE*(2. ODO+CE )/( (1.CD0+CE)**2) 

ENDIF 

Correction  of  the  cross  stream  wall  velocity  boundary 
condi t i on . 

40  SFICOR=SFIBC-SFI ( 1 ) 

Computation  of  the  outer  region  eddy  viscosity  according 
to  Clauser's  model. 

IF( PRTR. GT. 0 . ODO )  THEN 

EREF=0 . 0168D0*DSQRT( 2 .ODO*ZI(IS-1)) *DELSTR/VI SCL 
VISCT=VISCL* ( 1 . ODO +EREF* PRTR) 

ENDIF 

Turbulent  flow  strearawise  stretching  parameter. 
PRST=VISCL/VISCT 

Computation  of  the  coefficient  block  matrices. 

CALL  COEFFS ( AHAT, BHAT, CHAT, DHAT, S , P, Q ) 

CALL  KA(AHAT,BHAT,CHAT,A) 

CALL  MB (AHAT, BHAT, DHAT, 3) 

CALL  MC(C) 

CALL  MR(N,SFICCR,P,Q,S,R) 

Computation  of  the  matrices  used  in  the  Lower-Upper 
decomposition  of  the  block  tridiagonal  system. 

DO  50  1-1,3 
DO  45  J-1,3 

AL(I,J,1)=A(I,J,1) 

45  CONTINUE 
50  CONTINUE 

DO  100  K-2,NP1 
DO  75  I-l , 2 
DO  60  L=l,3 
DO  55  J=l,3 


onono  onn  nnon  nno 


DU(L, J)-AL( J,L,K-1) 
CONTINUE 
CONTINUE 
DO  65  J=l,3 

V( I , J,K) 
CONTINUE 
MS-1 
NS-3 
I  A- 3 
IDGT=0 

IMSL  matrix  solver. 


CALL  LEQT1F( DU,MS ,NS , lA, V, IDGT, WKAREA, lER) 

DO  70  J=l,3 

BET( I , J, K)=V( J, 1 ) 

70  CONTINUE 

75  CONTINUE 

DO  80  J=l,3 

BET( 3, J,K)=0.0D0 
80  CONTINUE 

DO  90  1=1,2 
DO  85  J=l,3 

AL(I, J,K)=A(I,J,K)-BET(I,3,K)*C( 3, J,K-1 ) 

85  CONTINUE 

90  CONTINUE 

DO  95  J=l,3 

AL(3,J,K)=A(3,J,K) 

95  CONTINUE 
100  CONTINUE 

The  backward  and  forward  substitutions  of  the  Lower 
Upper  decomposition  of  the  tridiagonal  system. 

CALL  FORW(N,BET,R,Y) 

CALL  BACK(N,Y,C,AL,Z) 

Update  of  the  station  computations. 

DO  110  K=1,NP1 

SFI(K)=SFI(K)+Z( 1,K) 

BFI(K)=BFI(K)+Z(2,K) 

GI{K)-GI(K)+Z( 3,K) 

110  CONTINUE 

Station  convergence  test  for  the  momentum  equation 
based  on  both  the  wall  velocity  gradient  correction 
and  the  outer  edge  velocity  gradient. 

IF(DABS(Z( 3,1) ) .GT.EPS)  THEN 
ITER( IS)-ITER( IS)+1 
GO  TO  40 
ELSE 

IF{DABS(GI(NP1 )) .GT.EPS)  THEN 


nnnn  non  nnnno  nnon  nnn  non  nnn 


Boundary  layer  extension  based  on  the  zero  edge 
velocity  gradient  condition. 

CALL  EXTENDI  ETA) 

NP1=N+1 

Update  of  the  normed  cross  stream  coordinate. 

DO  120  K=2,NP1 

YDEL( K)-ETA(K)/ETA(NP1 ) 

CONTINUE 

ITER( IS)=ITER(IS)+1 
GO  TO  40 
ENDIF 

Computation  of  the  void  fraction. 

IF( ( IS .GE. ISINJI ) .AND. ( ISINJF . NE . 0 ) )  THEN 
CALL  VOLFRC(XI ,XI1 , ZI ) 

Computation  of  the  buoyancy  force  field  according 
to  Archimedes  principle. 

CALL  ARCHIE(XI ,BUOY) 

Convergence  of  the  iteration  procedure  of  the  vhole 
system  based  on  the  convergence  of  the  wall  velocity 
gradient. 

SYSCON-DABS ( SYSCON-GI ( 1 ) ) 

IFISYSCON.GT.EPS)  THEN 
SYSCON-DABS (GI(1) ) 

ITER( IS ) -ITER( IS )+l 
GO  TO  40 
ENDIF 
ENDIF 

SYSCON-0 . ODO 

Computation  of  the  skin  friction. 

IF(IS.NE.l)  THEN 

CF(IS)-VISCL*GI(1) *DSQRT( 2 . ODO/ZI ( IS ) ) 

ENDIF 

Station  advance  and  specification  of  the  corresponding 
starting  solution. 

IS-IS+1 
ITER( IS )-l 
DO  130  K-1,NP1 
SFI2(K)-SFI1(K) 

SFIl ( K )-SFI ( K  ) 

BFI1(K)-BFI( K) 

GI1(K)-GI{K) 

XI2(K)-XI1(K) 


non 


XIl(K)=XI(K) 

SUM2(K)=*SUM1(K) 

SUM1(K)=SUM(K) 

BU0Y1{K)=BU0Y(K) 

130  CONTINUE 

Computation  of  the  boundary  layer  mass  flowrate. 

IF( ( IS .GT. ISINJI ) .AND. ( ISINJF.NE. 0 ) )  THEN 
PRQBL=0 . ODO 
DO  135  K=2,NP1 

CK=(1.0D0-XI(K) )*BFI(K) 

CK1=( 1 . ODO-XI ( K-1 ) ) *BFI ( K-1 ) 

PRQBL=PRQBL+( CK+CKl ) *DETA( K ) 

135  CONTINUE 

QBL( IS-1 )=PRQBL*DSQRT(ZI ( IS-1 )/2 . ODO ) 

ENDIF 
C 

C  The  end  of  the  computations  is  commanded  at  30  boundary 

C  layer  thicknesses  downstream  from  injection. 

C 

IF( ( (IS-l) .GT.ISINJF) .AND. ( I S INJ F . NE . 0 ) )  THEN 
POSEND=POSEND+DZI/(UE*ROUL*VISCTF) 

ENDIF 

IF(XDEL(IS-2) .GT.30.0D0)  GO  TO  150 
C 

C  Update  of  the  stceamwise  stepsize  during  transition. 

C 

IF( ( PRTR.GT. 0 . ODO ) .AND. ( PRTR . LT . 1 . ODO ) )  THEN 
DZIT=DZIC* ( 1 . ODO-O . 5*EREF*PRST* PRTR ) 

ELSE 

IF(PRTR.EQ.O.ODO)  THEN 
DZI=DZIC 
ELSE 

DZI=DZIC* ( 1 . ODO+0 . 8*EREF) 

ENDIF 

ENDIF 

C 

C  Update  of  the  streamwise  stepsize  during  injection. 

C 

IF( ( ISINJI .NE . 0 ) .AND. ( I S I N J F . NE . 0 ) )  THEN 
DZI=DZI/4 . ODO 
ENDIF 

IF( IS.GE. ( ISINJF+20) )  THEN 
DZI-4*DZI 
ENDIF 

ZI ( IS )-ZI ( IS-1 ) +DZI 

ALFA( IS ) -{ ZI (  IS  )+2: (  IS-1  ;  )/'ZZ  IS  :  -  Z  I  ;  IS-l  ' 

C 

C  Printing  statements. 

C 

IF(  (  (  IS-1  )  .GE.  ISINJF  )  .AND.  IS  IN’J.-  .N’E  .  :  .  THEN 

LP-IS-ISINJF-1 

IF(MOD{  DFLOATC  L?  )  ,  15  .  OdO  .  E  J.  I  .  IDC  ^  THEN 
L1-31+LP/15 


L2-51+LP/15 
L3=71+LP/15 
L4-91+LP/15 
WRITE( 21 , * ) 


XDEL( IS-2  ) 


Computation  of  the  velocity  profile  in  inner  variables. 

DO  140  K=2,NP1 

UPLUS( K)=BFI ( K)/USTAR 
YTRANS=DLOG( YPLUS( K) )/STRFAC 
WRITE (L4,*)  YTRANS,UPLUS( K) 

CONTINUE 
DO  145  K=1,NP1 
VEL=UE*BFI ( K ) 

VELGRD=ROUL*(UE**2 ) *GI ( K)/DSQRT( 2*ZI(IS-1)) 
WRITE(L1,*)  VEL,YDEL(K) 

WRITE(L2,*)  VELGRD, YDEL( K) 

WRITE(L3,*)  XI(K) ,YDEL(K) 

CONTINUE 
ENDIF 
ENDIF 
GO  TO  15 
ISEND=IS-3 
DO  155  K=2,ISEND 

REXLOG=DLOG(REX(K) )/STRFAC 

Computation  of  the  boundary  layer  thickness  and  momentum 
thickness  Reynolds  numbers. 

FRICK=1000*CF( K) 

RBL=ROUL*UE*THBL( K)/VISCL 
RMO=ROUL*UE*THMO( K )/VISCL 
RBLLOG=DLOG( REL)/STRFAC 
RMOLOG=DLOG ( RMO ) /STRFAC 
WRITE(120,*)  REXLOG, FRICK 
WRITE(121,*)  REXLOG, RBLLOG 
WRITE(122,*)  REXLOG, RMOLOG 
WRITE{123,*)  REXLOG, SHAPE(K) 

IF( K.GE. ISINJF)  THEN 
QRATIO=ROUG*QG/QBL( K) 

WRITE(124,*)  XDEL(K) ,QRATIO 
ENDIF 
CONTINUE 
ENDIF 


★  ★ 

*  SUBROUTINES  * 

it  ★ 


c 

C  Subroutine  FORW  performs  the  forward  substitution  of 

C  the  Lower-Upper  decomposition  of  the  block  tridiagonal 

C  system. 

C 


SUBROUTINE  FORW( N , BET , R , Y ) 

IMPLICIT  REAL*8  (A-H,0-Z) 

DIMENSION  Y(3,101),R(3,101) , BET ( 3, 3,101 ) 
NP1=N+1 
DO  10  J-1,3 
Y( J,l)-R( J,l) 

10  CONTINUE 

DO  40  K-2,NP1 
DO  30  I=«l,2 
Bl-0 . ODO 
DO  20  J=l,3 

B1=B1+BET( I, J,K)*Y( J,K-1) 

20  CONTINUE 

Y( I , K )-R( I , K )-Bl 
30  CONTINUE 

Y(3,K)-R(3,K) 

40  CONTINUE 
RETURN 
END 


C 

C  Subroutine  BACK  performs  the  backward  substitution  of 

C  the  Lower-Upper  decomposition  of  the  block  tridiagonal 

C  system. 

C 


SUBROUTINE  BACK ( N , Y , C , AL , Z ) 


Lli 


IMPLICIT  REAL*8  (A-H,0-Z) 

• 

DIMENSION  AL{3,3,101),C(3,3,101),Z(3,101),Y(3,101) 

DIMENSION  WKAREA( 10),V(3,1),DU(3,3) 

m 

.••Jv 

NPl-N+1 

DO  5  1=1,3 

V(I,1)=Y(I,NP1) 

5 

CONTINUE 

• 

DO  15  1=1,3 

J 

DO  10  J=l,3 

DU( I , J )=AL( I , J , NPl ) 

10 

CONTINUE 

-'A 

15 

CONTINUE 

'1^ 

MS-1 

NS  =  3 

# 

IA-3 

A: 

IDGT-0 

A 

c 

c 

IMSL  matrix  solver. 

c 

CALL  LEQT1F(DU,MS,NS,IA,V,IDGT,WKAREA, lER) 

• 

DO  20  1=1,3 

Z(I,NP1)=V(I,1) 

•A 

20 

CONTINUE 

DO  50  K-1,N 

»V„v 

Bl-0 . ODO 

•A 

DO  25  J=2 , 3 

• 

B1-B1+C( 3 , J,NP1-K) *Z( J,NP1+1-K) 

,N  ■ 

25 

CONTINUE 

Y( 3,NP1-K)=Y( 3 ,NP1-K)-B1 

DO  30  1=1,3 

V( I  ,1 )-Y( I ,NP1-K) 

30 

CONTINUE 

• 

DO  40  1=1,3 

DO  35  J=l,3 

DU( I , J )»AL( I , J,NP1-K) 

-  A 

35 

CONTINUE 

40 

CONTINUE 

A 

MS-1 

NS-3 

• 

IA-3 

IDGT-0 

. *•  ,*» 

c 

c 

IMSL  matrix  solver. 

c 

CALL  LEQT1F(DU,MS,NS,IA,V,IDGT,WKAREA, lER) 

DO  45  1=1,3 

Z{ I ,NP1-K )-V( 1,1) 

45 

CONTINUE 

» * '  .*» 

50 

CONTINUE 

v'‘.' 

RETURN 

END 

• 

J* 

I 


WF  IL*'.l^  IL"  IL*  H'W.UA'  LVJ. 


C 

C  Subroutine  COEFFS  computes  the  entries  of  the  matrices 

C  [A],  [B],  [C],  and  [W]  of  the  block  tridiagonal  system. 

C  The  system  arises  from  the  discretization  of  the 

C  parabolic  momentum  equation  with  Keller's  box  scheme. 

C 


SUBROUTINE  COEFFS ( AHAT , BHAT , CHAT , DHAT , S ,P,Q) 

IMPLICIT  REAL*8  (A-H,0-Z) 

DIMENSION  AHAT ( 101 ) , BHAT{ 101 ) ,DHAT( 101),GI(101),BFI{101) 
DIMENSION  P(101),Q(101),T(101),S(101) ,CHAT( 101 ) 

DIMENSION  SFI 1(101), BFI 1(101), GI 1(101), SFI2(101) 

DIMENSION  XI (101), XI 1(101) ,PXL( 101 ) , PXLl ( 101 ) ,BUOY( 101 ) 

DIMENSION  SFI ( 101 ) , BUOYl ( 101 ) ,DETA( 101), ALFA (1000) 

COMMON/BLKl/DETA , CE , N 

COMMON/BLK2/ALFA, IS 

COMMON/BLK3/SFI , BFI ,GI 

COMMON/BLK4/SFI1 ,BFI1 ,GIl 

COMMON/BLK5/PRTR, PRTRl , PRST,  PRSTl 

COMMON/BLK6/XI ,XI1 

COMMON/BLK8/PXL , PXLl 

COMMON/BLK14/BUOY, BUOYl 

NPl-N+1 

AI-ALFA( IS ) 

C 

C  Computation  of  the  mixing  length  based  on  the 

C  latest  station  update. 

C 

IF( PRTR .GT . 0 . ODO )  THEN 
CALL  MIXLEN( IREF, PRTR) 

ENDIF 

DO  10  K=2,NP1 
C 

C  Influence  coefficient  of  the  inner  eddy  viscosity. 

C 

IF( PRTR. GT. 0 . ODO )  THEN 
IF( ( K-1 ) .LE. IREF)  THEN 
PINFl-2 . ODO 
ELSE 

PINFl-1 . ODO 
ENDIF 

IF( K . LE . IREF )  THEN 
PINF-2 . ODO 
ELSE 

PINF-1 . ODO 
ENDIF 
ENDIF 

DE-DETA( K ) 

ASF-(SFI(K)+SFI(K-1))/2.0D0 
ASF1-(SFI1(K)+SFI1(X-1) )/2.0D0 
ABF-(BFI(K)+BFI(K-1) )/2.0D0 
ABF1-(BFI1(K)+BFI1(K-1))/2.0D0 


o  o  o 


AG-(GI(K)+GI (K-1 ) )/2.0D0 
AG1=(GI1(K)+GI1(K-1) )/2.0D0 
AX-(XI ( K)+XI ( K-1 ) )/2.0d0 
AX1=(XI1(K)+XI1(K-1) )/2.0D0 
AXD-XI(K)-XI(K-1) 

AX1D=XI1(K)-XI1(K-1) 

PVF=1 . 0D0-( AX+AXl )/2 . ODO 
PVK=1 .0D0+(0.9+2.6*XI(K))*XI(K) 

PVK1=1 . ODO+( 0.9+2.6*XI(K-l))*XI(K-l) 

PV1K=1 . 0D0+( 0.9+2.6*XIl(K))*XIl(K) 

PV1K1=«1 .0D0+(0.9  +  2.6*XIl(K-l))*XIl(K-l) 

PLKG-1 . 0D0+( PINF/{ l.ODO-XI(K) ) )*PXL(K) 

PLKG1  =  1 .0D0+(PINFl/(l. ODO-XI ( K-1  )  )  ) *PXL ( K-1  ) 
PLK-1.0D0+(1.0d0/(1.0D0-XI(K) ) )*PXL(K) 

PLKl-1 . 0D0+( 1 . 0D0/( 1 . ODO-XI(K-l ) ) ) *PXL( K-1 ) 

PL1K=«1  .  0D0+(  1 . 0D0/(  1 .  ODO-XIK  K)  )  )  *PXL1  {  K) 
PLlKl-1.0D0+(l. 0D0/( 1 . ODO-XIK  K-1 )  )  ) *PXLl ( K-1  ) 

P( K)--SFI ( K )+DE*ABF  +  SFI ( K-1  ) 

Q(K)=-BFI(K)+DE*AG+BFI ( K-1  ) 

AHAT( K ) -DE*PVF*AI *ABF 

BHAT(  K)='-DE*PVF*  (AG+AI*  ( AG+AGl  >  )/2  .  ODO 

CHAT( K)=-PVF*DE* ( ASF-AI* (ASFl-ASF ) )/2 . ODO-PRST* 

>(AX+( 1 . ODO-AX) *PLKG) *PVK 

DHAT( K)=— PVF*DE* ( ASF-AI*(ASF1-ASF)  )/2 , ODO  +  PRSTl* 
>(AXl+( 1 . ODO-AXl ) *PLKG1 ) *PVK1 

51- DE*(AI*(ABF**2-A3F1**2+AG*ASF1-ASF*AG1 )-( 1 . ODO+AI ) * 
>ASF*AG-( 1 . ODO-AI ) *AG1*ASF1 ) 

52- AX*PRST* ( PVK*GI ( K ) -?VK1*GI ( K-1 )  )  + 

>AX1*PRST1*( PV1K*GI1 ( K)-PV1K1*GI1 ( K-1 ) ) 

S3=(l . ODO-AX) * ( PVK*PLK*GI ( K ) -PVK 1  * PLKl *GI ( K-1 )  ) * 
>PRST+{ 1 . ODO-AXl ) *( PV1K*PL1K*GI1 ( K)-PV1K1*PL1K1*GI1 ( K-1 ) ) 
>  *PRST1 

Computation  of  the  buoyancy  generated  source  term 
of  the  momentum  equation. 

S4=DE* ( BUOY( K ) +BUOY( K-1 )+BUCYl ( K ) +BUOY1 (K-1) )/2.0D0 
S( K) — S1*PVF+S2+S3+S4 
10  CONTINUE 
RETURN 
END 


Subroutine  MA  computes  the  matrix  [A]  of  the  discretized 
momentum  equation. 


SUBROUTINE  NA ( AHAT , BHAT , CHAT , A ) 

IMPLICIT  REAL*8  (A-H,0-Z) 

DIMENSION  DETA( 101 ) , A( 3 , 3 , 101 ) 

DIMENSION  AHAT( 101  )  , CHAT ( 101 )  ,BHAT( 101 ) 

COMMON/BLKl/DETA,CE,N 

NPl-N+1 

A( 1 , 1 , 1 )-l . ODO 

A(2,2,1)-1.0D0 

A( 3 , 2 , 1 )=-l . ODO 

A( 3 , 3 , 1 )=-DETA( 2 )/2 . ODO 

A(1,2,1)=O.ODO 

A(  2 , 1 , 1  )=0 . ODO 

A(  2 , 3 , 1  )=>0  .  ODO 

A( 1 , 3 , 1  )=0 . ODO 

A(  3,1,1)«0.0D0 

DO  10  K-2,N 

A( 1 , 1 , K )=1 . ODO 

A( 1 , 2 , K )--DETA( K )/2 . ODO 

A(2,1,K)=BHAT(K) 

A(2,2,K)-AHAT(K) 

A(  2 , 3 , K)-CHAT( K ) 

A( 3 , 1 , K )-0 . ODO 
A(3,2,K)— l.ODO 
A( 3 , 3 ,K)=— DETA( K+1 )/2 .ODO 
A(  1 ,  3  ,  K)=«0.  ODO 
10  CONTINUE 

A( 1 , 1 ,NP1 )=1 . ODO 

A( 1 , 2 ,NP1 )--DETA( NPl )/2 . ODO 

A( 1 , 3 ,NP1 )=0 . ODO 

A( 2, 1 ,NPl )-BHAT(NPl ) 

A(2,2,NP1)=AHAT(NP1) 

A(2,3,NP1)=CHAT(NP1) 

A(  3 , 1 ,NP1 )-0 . ODO 
A( 3 , 2 ,NP1 ) =1 . ODO 
A( 3 , 3 ,NP1 )=0 .ODO 
RETURN 
END 


Subroutine  MB  computes  the  matrix  [B]  of  the  discretized 
momentum  equation. 


SUBROUTINE  MB ( AHAT , BHAT , DHAT , B ) 

IMPLICIT  REAL*8  (A-H,0-Z) 

DIMENSION  B( 3 , 3  ,  101  )  ,DETA( 101 ) 

DIMENSION  AHAT( 101 ) ,BHAT( 101 ) ,DHAT( 101 ) 


COMMON/BLKl/DETA , CE  ,  N 
NP1=N+1 
DO  10  K=2,NP1 
B(1,1,K)=-1.0D0 
B( 1 , 2,K)=-DETA(K)/2.0D0 
B  (  1 , 3 , K)=0 , ODO 
B(  2 , 1 , K)=BHAT( K) 
B(2,2,K)=AHAT(K) 

B( 2 , 3 , K)=DHAT( K) 

B( 3 , 1 , K)-0 . ODO 
B( 3 , 2 , K)=0  .  ODO 
B{ 3 , 3 , K)=0 . ODO 
10  CONTINUE 
RETURN 
END 


Subroutine  MC  computes  the  matrix  IC}  of  the  discretized 
momentum  equation. 


SUBROUTINE  MC ( C ) 

IMPLICIT  REAL*8  (A-H,0-Z) 
DIMENSION  C( 3 , 3 , 101 ) ,DETA( 101  ) 
COMMON/BLKl/DETA , CE , N 
DO  10  K=1,N 

C( 3,1,K)=O.ODO 
C( 3 , 2 , K)-l . ODO 
C( 3 , 3 , K)=-DETA(K+1 )/2 . ODO 
C( 1,1,K)=0.0D0 
C(1,2,K)=0.0D0 
C( 1 , 3 , K ) -0 . ODO 
C(2,1,K)-0.0D0 
C(2,2,K)=0.0D0 
C(2,3,K)-0.0D0 
10  CONTINUE 
RETURN 
END 


Subroutine  MR  computes  the  vector  [W],  the  source  term 
of  the  discretized  momentum  equation. 


U  U  u 


SUBROUTINE  MR ( N , S F I  COR , P , Q , S , R ) 

IMPLICIT  REAL*8  (A-H,0-Z) 

DIMENSION  R(3,101),P(101),Q(101),S(101) 

NP1=N+1 

R( 1 , 1 )=SFICOR 

R( 2 , 1 )=0 . ODO 

R(3,1)=Q(2) 

DO  10  K=2,N 
R( 1,K)=P(K) 

R( 2,K)=S(K) 

R( 3,K)=Q(K+1 ) 

10  CONTINUE 

R( 1 ,NP1 )=P(NP1 ) 

R( 2 ,NP1 )=S(NP1 ) 

R( 3,NP1)=0.0D0 

RETURN 

END 


Subroutine  EXTEND  performs  the  extension  of  th 
computational  domain  when  necessary. 


SUBROUTINE  EXTEND(ETA) 

IMPLICIT  REAL*8(A-H,0-Z) 

DIMENSION  DETA( 101),SFI(101),BFI(101),GI(101) 

DIMENSION  SFI1(101),BFI1(101),XI(101),XI1(101) 

DIMENSION  PXL( 101  )  , PXLl ( 101 ) , ETA( 101 )  , SUM{ 101  ) 

DIMENSION  XI2 ( 101  )  , BUOY( 101 ) ,BUOYl ( 101 ) 

DIMENSION  GIl ( 101 ) , SUM2 ( 101 ) , SUMl ( 101 ) 

COMMON/BLKl/DETA , CE , N 

COMMON/BLK3/SFI ,BFI ,GI 

COMMON/BLK4/SFI1 ,BFI1 ,GI1 

COMMON/BLK6/XI ,XIl 

COMMON/BLK8/PXL , PXLl 

COMMON/BLKll/SUM , SUMl , SUM2 

COMMON/BLK13/SFI2 , XI2 

COMMON/BLK14/BUOY, BUOYl 

NREF-N 

N-N+2 


DO  10  J=N,NP1 

DETA( J)=DETA( J-1)*CE 
ETA( J)=ETA( J-1)+DETA( J) 
PXL( J)=PXL( J-1) 

PXLl( J)=PXL1( J-1) 

XI( J)=XI( J-1) 

XIl( J)=XIl( J-1) 

XI2( J)=XI2( J-1) 

SUM( J)=0.0D0 
SUMK  J)=O.ODO 
SUM2 ( J ) =0 . OdO 
BUOY( J ) =0 . ODO 


BUOYl ( J )=0 . ODO 
CONTINUE 
DO  20  J=NMl,NPl 

SLOPE=( ETA( J)-ETA(NREF) )/{ETA(NPl )-ETA(NREF) ) 

BFIK J)=BFIl(NREF)+( 1 . ODO-BFI 1 ( NREF ) )* SLOPE 
SFI1(J)=SFI1( J-1 )+DETA( J) * (BFI1(J)+BFI1(J-1))/2.0D0 
BFI ( J )=BFI (NREF)  +  ( 1 . ODO-BFI ( NREF )  ) *  SLOPE 
SFI ( J  )=SFI ( J-1 )+DETA( J ) * ( BFI ( J ) +BFI ( J-1 )  ) /2 . ODO 
SFI2( J)=SFI2( J-1)+DETA( J)*(BFI{ J)+BFI( J-1) )/2.0D0 
GI{ J)=0.0D0 
GIK  J)=0.0D0 
CONTINUE 
RETURN 
END 


m 


Subroutine  MIXLEN  computes  the  wall  function  of 
the  turbulence  model  through  a  modified  Van  Driest 
type  expression. 


SUBROUTINE  MIXLEN ( I  REF , PRTR ) 

IMPLICIT  REAL*8(A-H,0-Z) 

DIMENSION  PXL( 101 ) , PXLl (101),SFI(101),BFI(101) 
DIMENSION  GI { 101 ) ,ETA( 101),XI{101),XIl(101),ZI{1000 
DIMENSION  SFI1(101),BFI1(101),GI1(101), ALFA (1000) 
DIMENSION  YPLUS(lOl) ,DETA(101) 

COMMON/BLKl/DETA , CE , N 
COMMON/BLK2/ALFA, IS 
COMMON/BLK3/SFI ,BFI ,GI 
COMMON/BLK4/SFIl,BFIl ,GI1 
COMMON/BLK6/XI ,XI1 
COMMON/BLK7/ZI , ETA 
COMMON/BLK8/PXL , PXLl 

COMMON/BLK9/ISINJI , EREF , VI SCL , VI SCT 


C0MM0N/BLK12/YPLUS , USTAR 
NP1=N+1 

DZI=ZI ( IS )-ZI ( IS-1 ) 

DO  5  K-2,NP1 

PXLl ( K )=PXL( K ) 

5  CONTINUE 

PRZI=DSQRT( 2*ZI ( IS ) ) 

USTAR=DSQRT( DABS ( GI ( 1 ) ) *VISCL/PRZI ) 

VWALL=-VISCT* (SFI(1)+2*ZI(IS)*(SFI(1)-SFI1(1) )/DZI )/PRZI 

Wall  function  injection  modification. 

IF( IS .LE. ISINJI )  THEN 
ZETA-0 . 5D0 
ELSE 

ZETA=1 , 5D0 
ENDIF 

APLUS  =  26 . 0D0/{  1 .  ODO  +  5 . 9*\nvALL/USTAR  ) 

DO  20  K=2,NP1 

RA1=(1.0D0-XI(K) )/( 1 . ODO-XI ( 1 ) ) 

RA2=( 1 . 0D0+( 0 . 9+2 . 6*XI ( K) )*XI(K) )/ 
>(1.0D0+(0.9+2.6*XI(l))*XI(l)) 

RA3=(1. ODO-XI (K) )/(1.0DO+(0.9+2.6*XI(K) )*XI(K) ) 
YPLUS(K)=ETA(K) *DSQRT(PRZI*RA1*RA3*DABS(GI( 1) )/ 

X VISCL*RA2 ) ) 

P0WER=-YPLUS1K)* I ( ( RA2/RA1 ) *DABS ( GI ( K ) /GI { 1 ) ) )** 
>ZETA)/APLUS 

EII  =  0 . 41*ETA( K) *{ 1 .0D0-DEXP( POWER)  ) 
PXL(K)=PRZI*DABS(GI(K)  ) * ( El I  *  *  2 ) /VI SCL 
20  CONTINUE 
E0I=EREF 
DO  25  K=2,NP1 

IF(PXL(K) .LT.EOI )  THEN 
IREF=K 
ELSE 

GO  TO  30 
ENDIF 

25  CONTINUE 
30  IREFP1»IREF+1 

DO  40  K-IREFP1,NP1 
PXL( K  )  =E0I 
40  CONTINUE 

Modification  of  the  mixing  length  during 
transition . 

IF( PRTR. LT. 1 . ODO )  THEN 
DO  50  K-2,NP1 

PXL( K)-PRTR*PXL( K ) 

50  CONTINUE 
ENDIF 
RETURN 


n  n  n  n 


Subroutine  BLPR  computes  boundary  layer  integral 
parameters,  the  boundary  layer  displacement 
thickness  and  the  momentum  thickness. 


SUBROUTINE  BLPR ( DELSTR , TETSTR ) 

IMPLICIT  REAL*8(A-H,0-Z) 

DIMENSION  ETA( 101 ) ,DETA( 101), 21(1000), XI{101) 

DIMENSION  SFI(101),BFI{101),GI(101),XI1(101) 

COMMON/BL.Kl/DETA,  CE,N 

COMMON/BLK3/SFI ,BFI ,GI 

COMMON/BLK6/XI ,XI1 

NP1=N+1 

DELSTR=0 . ODO 

TETSTR=0 . ODO 

DO  10  K=2,NP1 

AKl-( 1  .  ODO-XI ( K-1 )  ) * ( 1 . ODO-BFI ( K-1  )  ) 

AK-( 1 . ODO-XI(K)  )  *{ 1 . ODO-BFI ( K) ) 

BK1=AK1*BFI ( K-1 ) 

BK=AK*BFI ( K ) 

DELSTR=DELSTR+ ( AKl+AK ) *DETA( K ) 
TETSTR=TETSTR+ ( BKl+BK ) *DETA( K ) 

CONTINUE 

DELSTR=DELSTR/2 . ODO 
TETSTR=TETSTR/2 . OdO 
RETURN 
END 


Subroutine  VOLFRC  computes  computes  the  volume  fractrcn 
according  to  a  hyperbolic  transport  equation. 


SUBROUTINE  VOLFRC ( X I , XI 1 , 2 1 
IMPLICIT  REAL*8 { A-H , 0-2 ) 
DIMENSION  DETA{ 101 ) , XI ( 101 ) 
DIMENSION  SFI ( 101 ) ,BFI ( 101 ) 
DIMENSION  SFI2( 101 ) ,XI2(101 
COMMON/BLKl/DETA,CE,N 
COMMON/BLK2/ALFA, IS 
COMMON/BLK3/SFI ,BFI ,GI 
COMMON/BLK4/SFI1 ,BFI1 ,GI1 
C0MM0N/BLK13/SFI2,XI2 


XIKlOl),  ALFA  (1000), GIKlOl) 
GI(101),SFI1(101),BFI1(101) 
,21(1000) 


EPS=1 . OD-5 
NPl=N+l 
BI=ZI ( IS ) 

DZI  =  ZI ( IS )-ZI ( IS-1  ) 

DZIl=ZI(IS-l)-ZI(IS-2) 

DE=DZI/DZI1 

Computation  of  the  first  node  of  the  void 
fraction  profile  after  injection. 

IF( IS.GT.ISINJF)  THEN 

PRI=( 2*DE+1 . ODO )/( 1 . ODO+DE) 

PRIl=-( 1 . ODO+DE) 

PRI2=( DE**2 )/( 1 . ODO+DE) 

DIR=2*BI*DETA( 2 ) *BFI ( 2 ) *(PRIl*XIl ( 2 )+PRl2*Xl2  (  2 )  ) 
COEFl  =  PRI*SFI ( 2 )+PRIl*SFIl ( 2 ) +PRI 2  * SFI 2 ( 2 ) 
ESQ=2*BI*DETA( 2 ) *BFI ( 2 ) *PRI-DZI *SFI ( 2 ) -2 *BI *COEFl 
XI(2)=DIR/ESQ 
END  IF 

Computation  of  the  void  fraction  profile. 

5  DO  10  K=3,NP1 

BIK=2*DETA( K) *BI*BFI (K) 

CIK=DZI*(  1 . ODO  +  DE) *SFI (K)+2*BI*( SFI2 ( K)  *DE**2- 
>SFI1(K)*(1. ODO+DE) **2+( 2*DE+1 . 0D0)*SFI(K) ) 

COEFl=CIK/( ( 1 . 0D0+2*DE)*BIK-CIK) 

COEF2=BIK/( ( 1 .0D0+2*DE)*BIK-CIK) 

XI  {  K)=-COEFl*XI ( K-1 )-C0EF2* 

>(XI2( K) *DE**2-XI1( K) *( 1 .0D0  +  DE)**2  ) 

10  CONTINUE 

Convergence  procedure  of  the  void  fraction  profile 
based  on  the  zero  wall  void  fraction  gradient. 

IF{ IS . LE. ISINJF )  THEN 

PRXIl=-(2.ODO+CE)/(DETA(2)*(CE+1.0DO) ) 

PRXI2=( 1 . 0D0+CE)/DETA{ 3 ) 

PRXI3=-1 . 0D0/( DETA( 3 ) * ( 1 . ODO+CE ) ) 

XIWGRD=XI ( 1 ) *PRXI1+XI ( 2 ) *PRXI2+XI ( 3 ) *PRXI3 
IF{DABS(XIWGRD) .GT.EPS)  THEN 

XI(2)  =  (XI(1)*CE*(2.0DO+CE)+XI(3) )/(  ( 1 . ODO  +  CE ) *  *  2 ) 
GO  TO  5 
ENDIF 
ENDIF 
RETURN 
END 


Subroutine  ARCHIE  computes  the  source  term  of  the 


momentum  equation  due  to  the  buoyancy  force  field 
according  to  Archimedes  principle. 


SUBROUTINE  ARCHI £ ( XI , BUOY ) 

IMPLICIT  REAL*8(A-H,0-Z) 

DIMENSION  SUM( 101 )  , SUMl ( 101  )  , DETA( 1 0 1 )  , ALr A (  1  0 0 0  ) 
DIMENSION  ZI  ( 1000  )  , ETA( 101),XI(101), BUOY( 101 ) 
DIMENSION  XII ( 101 )  , SUM2 ( 101  ) 

COMMON/BLKl/DETA , CE , N 
COMMON/BLK2/ALFA, IS 
COMMON/BLK7/ZI , ETA 
COMMON/BLKIO/GRAV, ROUL , UE 
COMMON/BLKl 1/SUM , SUMl , SUM2 
NP1=N+1 

DZI-ZI(IS)-ZI(IS-1) 

DZIl-ZI(IS-l)-ZI(IS-2) 

DE-DZI/DZIl 

RAX-2*ZI ( IS )/DZI 

FAC1*( 2*DE+1 . ODO )/( 1 .ODO+DE) 

FAC2='-(  1  ,  ODO  +  DE) 

FAC3-(DE**2)/( 1. ODO+DE) 

FAC4-( 2*CE+1 . ODO )/( 1 . ODO+CE) 

FAC5=- ( 1 . ODO+CE ) 

FAC6-( CE**2 )/{ 1 . ODO+CE) 

DO  10  K-1,N 

SUM(N+l-K)=SUM(N+2-K)+DETA(N+2-K) *(XI(N+1-K)+ 

>XI (N+2-K) )/2 . ODO 
10  CONTINUE 

DISTX=DSQRT( 2*ZI { IS ) ) 

FACTOR-DISTX*GRAV/( ROUL* { UE**3 ) ) 

BUOY( 1 ) -FACTOR* ( SUM ( 1 )+RAX*( SUM( 1 ) *FACl+SUMl ( 1 ) *FAC2- 
>SUM2{ 1 ) *FAC3 ) ) 

TERM1-SUM( 2 )-SUM( 1 ) 

TERM2=SUM( 2 ) *FACl+SUMl( 2) *FAC2+SUM2( 2) *FAC3 
BUOY( 2 ) -FACTOR* ( SUMl 2 ) +RAX*TERM2-TERMl ) 

DO  20  K-3,NP1 

RAY-ETA( K )/DETA( K ) 

TERM1-SUM( K) *  FAC 4  + SUMl K-1 ) *  FAC 5  + SUM  I K-2 ) *FAC6 
TERM2-SUMI K) *FAC1+SUM1 [ K ) *AFC2+SUM2 | K ) *FAC3 
BUOY! K ) -SUMl K ) +RAX*TERM2-RAY*TERMl 
BUOY! K ) -BUOY! K) *FACTOR 
20  CONTINUE 
RETURN 
END 
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